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SECTION  1 


INTRODUCTION  AND  SUMMARY 

During  this  contract,  three  papers  were  published  and  two  others 
are  in  print.  A sixth  paper  is  now  being  prepared  for  publication. 

These  papers  (listed  in  chronological  order  as  Refs.  1 through  6)  are 
included  as  Appendices  1 through  6. 

These  papers  fit  into  two  categories:  (1)  the  study  of  the  scalar- 

1 3 

product  (SP)  expression  of  the  boundary  free  energy  ’ and  (2)  the  study 

of  the  structure  of  the  antiphase  (APB)  and  intemhase  (IPB)  boundaries 

2 A 5 6 

and  grain  boundaries  using  the  sum  method  ’ ’ ’ rather  than  using  the 
SP  method. 

If  used  with  care,  the  SP  method  is  simpler  than  the  sum  method  in 
obtaining  the  boundary  free  energy  value.  Appendix  1 concludes  that  the 
SP  method  is  reliable,  and  Appendix  3 extends  the  original  SP  formulation 
to  the  case  of  long-range  interaction  energy.  However,  during  the  three- 
year  contract  period,  the  SP  method  has  not  been  used  to  its  full  capacity. 

The  sum  method  has  the  advantage  that  it  gives  the  profile  (which 
the  SP  method  cannot  give)  across  the  boundary;  however,  the  computation 
is  more  lengthy  than  with  the  SP  method.  Appendix  2 introduces  a crucial 
concept  that  opens  up  a way  of  using  the  sum  method  to  the  boundary 
structure,  and  the  procedure  is  shown  with  an  example  in  a b.c.c.  ordered 
structure. 

The  method  developed  in  Appendix  2 is  put  to  full  use  in  Appendices  4 
and  5,  in  which  the  APBs  within  the  ordered  Cu^Au  phase  and  the  IPBs 
between  the  Cu^Au  phase  and  the  disordered  phase  are  studied.  This  work 
is  a revival  of  cooperation  on  a similar  subject  with  Dr.  John  W.  Cahn 
(now  at  the  National  Bureau  of  Standards).  Appendices  4,  5,  and  6 rely 
heavily  on  Dr.  Cahn's  knowledge  concerning  the  metallurgical  aspects  of 
the  problems  considered. 

Three  new  findings  came  out  of  the  study  of  the  Cu^Au  phase  boundar- 
ies. The  first  is  the  behavior  of  the  APB  free  energy  a,  which  shows 
a maximum  at  a temperature  of  about  0.6  of  the  disordering  temperature 


i 


is  the  discovery  of  a series  of  second-order  phase  transitions  within 
the  APB  slightly  below  T . 

The  structure  of  grain  boundaries  (GBs)  was  studied  during  the 
contract.  Although  work  in  this  area  is  not  yet  complete,  preliminary 
findings  are  being  written  up  for  publication  (probably  in  Physical 
Review  Letters).  The  main  discoveries  are  first  that  distinct  dis- 
tinguishable low-  and  high-temperature  structures  of  the  GBs  exist,  and 
second  that  the  GB  is  completely  wet  at  the  melting  point. 


SECTION  2 


THE  NATIONAL  ITERATION  METHOD  WITH  CONSTRAINTS 

The  natural  iteration  method  (NIM)  was  introduced  some  time  ago^  in 
solving  high-order  simultaneous  algebraic  equations  appearing  in  the 
cluster-variation  method  (CVM).  As  NIM  applications  increased,  it 
became  necessary  to  treat  cases  in  which  subsidiary  conditions  were 
imposed.  Appendix  1 shows  that  these  subsidiary  conditions  can  be  treated 
based  on  a concept  similar  to  the  original  NIM.  Within  each  iteration 
step  originally  designed  (which  we  call  the  "major"  iteration),  we  do 
what  we  call  "minor"  iterations  to  satisfy  subsidiary  conditions.  In 
Ref.  7 we  proved  that  the  major  iterations  always  converge;  in  contrast, 
we  have  not  proved  that  the  minor  iterations  converge.  However,  in  all 
the  cases  during  the  last  three  years  in  which  we  used  the  minor  itera- 
tions, they  never  failed  to  converge. 

When  subsidiary  conditions  exist,  we  can  classify  the  variables  as 
independent  or  dependent.  One  way  of  solving  the  equilibrium  state  is 
to  minimize  the  free  energy  with  respect  to  the  independent  variables 
(rather  than  using  Lagrange  multipliers  for  subsidiary  conditions). 

When  this  is  done,  the  resulting  simultaneous  equations  are  not  of  the 
form  for  which  the  NIM  is  applicable.  In  such  a case,  the  Newton- 
Raphson  (N-R)  method  can  certainly  be  applied,  but  it  is  often  quite 
time  consuming  to  find  the  right  initial  guess  of  variables  for  the  N-R 
method,  since  otherwise  the  N-R  method  does  not  converge  to  the  desired 
equilibrium  state.  In  the  boundary  studies  described  below  in 
Appendices  2,4,5,  and  6,  it  is  practically  impossible  to  use  the  N-R 
method  because  there  are  several  thousand  independent  variables. 
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SECTION  3 


THE  SP  FORMULATION  OF  THE  BOUNDARY  ENERGY 


The  scalar-product  (SP)  formulation  of  the  boundary  free  energy, 

g 

first  proposed  by  Clayton  and  Woodbury,  calculates  the  boundary  free 
energy  a using  the  formula 


exp(-Ao/RT) 


-£[P1<V>  P2(V)’ 

v 


1/2 


(1) 


! 


where  A is  the  sectional  area  parallel  to  the  boundary.  This  formula  is 

N 

written  for  a lattice  structure  having  phase  1 on  the  left,  phase  2 on 

the  right,  and  the  boundary  in  between. 

We  consider  a lattice  plane  parallel  to  the  boundary  but  far  away 

from  it,  and  well  inside  the  bulk  phase  i (i  = 1 or  2).  A configuration 

within  the  plane  is  denoted  by  V,  and  P^(v)  is  the  probability  of  finding 

the  configuration  v in  the  plane  inside  the  bulk  phase  i (i  = 1 or  2). 

Since  the  right  side  of  Eq.  1 has  the  form  of  the  SP  of  two  vectors 

1/2 

[P1(v)J  with  i - 1 and  2,  we  call  Eq.  1 the  SP  formula  of  the  boundary 
free  energy  a. 

■ '■ 

Eq.  1 is  noteworthy  in  that  0 can  be  calculated  by  knowing  only  the 
properties  of  the  two  bulk  phases.  Expression  1 is  reasonable  since, 
when  phases  1 and  2 are  identical,  a vanishes  because  P^(v)  is  normalized 
to  unity  and  since,  when  the  two  bulk  phases  are  very  different  in  their 
properties,  0 is  large. 

Eq.  1 is  rigorous  when  V is  for  the  configuration  of  an  infinitely 
wide  plane  parallel  to  the  boundary.  In  Eq.  1,  however,  approximations 
must  be  introduced.  In  Ref.  9,  we  had  calculated  Eq.  1 using  the  general 
CVM  approach.  Another  feature  of  Eq.  1 is  that  a rigorous  proof  was 

g 

lacking;  Clayton  and  Woodbury's  proof  was  not  sufficient,  and  our  attempt 
in  Ref.  9 still  lacked  mathematical  rigor.  Therefore,  to  use  Eq.  1 with 
enough  confidence,  we  felt  it  necessary  to  use  better  approximations  in 
the  CVM  and  to  compare  the  result  with  Onsager's^  rigorous  result  in  a 
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two-dimensional  square  lattice  Ising  model.  In  Ref.  9,  we  used  the  pair 
approximation  and  the  square-angle  approximateion  of  the  CVM.  In 
Appendix  1,  we  improved  the  treatment  using  the  double  square  r or  the 
(l.O)-direction  boundary  and  the  W-shaped  cluster  for  the  (1, Indirection 
boundary,  both  for  the  two-dimensional  square  Ising  lattice.  Since  the 
results  are  converging  to  Onsager's  rigorous  result,  we  concluded  in 
Appendix  1 that  we  can  use  the  general  expression  (Eq.  1)  with 
conf idence. 

The  general  idea  behind  the  derivation  of  Eq.  1 is  given  below. 

Consider  a three-dimensional  system  as  made  up  of  two-dimensionally 

large  planes  stacked  on  top  of  each  other.  Write  the  free  energy  F of 

the  entire  system  using  tin  CVM  in  terms  of  the  probability  variables 

for  two  adjoining  planes  where  n indicates  the  location  of  the 

planes.  Minimizing  F vith  respect  to  P (.U,v)  yields  a relation  that 

n 1/2 

expresses  P^Cp.v)  as  proportional  to  [ Pn (y ) pn+^(v)i  • If  the  system 
is  homogeneous  (i.e.,  no  boundary  in  it),  this  procedure  leads  to  the 
eigenvalue  formulation  with  which  Onsager  started. ^ The  details  of 
the  derivation  of  Eq.  1 from  this  point  on  are  presented  in  Refs.  3 and  9. 


SECTION  4 


THE  SP  EXPRESSION  FOR  LONG-RANGE  INTERACTION 


• Eq.  1 considers  one  lattice  plane  in  each  bulk  phase.  This 
expression  is  good  when  the  interaction  potential  is  of  the  nearest- 
neighbor  type.  But  it  must  be  modified  when  the  interaction  potential 
goes  beyond  the  nearest  neighbor.  This  problem  is  worked  out  in 
Appendix  3.  The  resultant  formula  modifies  Eq.  1 as 


exp(-Ao/RT)  » ^ ^Pi(S,i,V2»  * * * P2(V1’V2'  (2) 

x exp|a^(v^,\>2»  . . . ,v^)  - a2(vi,v2’  "**,vk^  ’ 

where  P . . . ,v  ) is  the  extension  of  P . (v)  and  is  the  probability 

that  consecutive  1,2,  . . . ,k  planes  inside  the  bulk  phase  i take  the 

configurations  *‘*,vk*  var^a^^e  a^^v^»v2  * * ' ,vk^  3 

Lagrange  multiplier  to  guarantee  continuity  of  the  form 


^ ^P.^  (p , , ' * ' >vfc)  * 

V 5 

Eq.  2 can  take  into  account  the  long-range  interactions  up  to  the 
interaction  between  the  lattice  1 and  the  lattice  (k  + 1).  Eq.  2 reduces 
to  Eq.  1 when  k - 1 and  when  the  symmetry  of  the  lattice 


P^U.V)  = P1(V,p) 


(4) 


holds. 


11 


The  general  expression,  Eq.  2,  was  tested  for  the  two-dimensional 
Islng  model  using  a 3 x 2 cluster  (i.e.,  a double-square  cluster  made  of 
six  lattice  points)  with  the  "3"-side  perpendicular  to  the  boundary.3 
The  result  agrees  well  wUh  that  of  the  double  square  in  Appendix  1 when 
the  a terms  in  Eq.  2 are  Included. 
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SECTION  5 


THE  SUM  METHOD  OF  CALCULATING  BOUNDARIES 

For  lack  of  a better  name,  we  call  this  method  the  sum  method  in 

contrast  to  the  SP  method.  Both  methods  are  based  on  the  CVM  mode  of 

thinking, but  differ  as  to  the  step  at  which  and  the  method  by  which 

they  introduce  the  approximation  into  the  formulation.  The  two  methods 

produce  similar  but  different  results.  Although  the  SP  method  has  the 

advantage  that  it  can  define  the  transition  point  within  the  boundary 
3 

more  clearly,  it  cannot  give  the  information  about  the  structure  across 
the  boundary.  The  sum  method  can  give  the  latter  information,  but 
usually  at  the  cost  of  more  computer  time.  Sometimes  the  approximation 
used  in  reducing  the  SP  formulation  to  the  tractable  level  makes  the 
results  unacceptable  on  a physical  basis  (such  a case  is  discui-sed  in 
Appendix  5) . A similar  trouble  has  not  been  encountered  with  the 
approximate  treatments  of  the  sum  method. 

12  13 

When  the  sum  method  was  used  more  than  ten  years  ago,  ’ the  NIM 
was  unknown.  After  the  NIM  had  been  devised  in  1974,^  we  tried  to 
apply  it  to  the  boundary  structure.  However,  a difficulty  was  encoun- 
tered in  handling  the  normalization  condition  for  individual  planes 

14 

(parallel  to  the  boundary)  until  the  work  of  Weeks  and  Gilmer  was 
noticed . 

In  Appendix  2,  we  treat  the  boundary  between  the  + spin  phase  and 

the  - spin  phase  (of  the  Ising  model)  in  a b.c.c.  structure  using  a 

tetrahedron  (irregular)  as  the  basic  cluster.  Appendix  2 presents  how 

14 

the  idea  of  Weeks  and  Gilmer  can  be  incorporated  in  the  CVM-based 
treatment.  The  advantages  of  the  latter  over  the  original  method  in 
Ref.  14  are  that  the  excess  free  energy  is  obtained  with  ease  from  the 
treatment  and  that  it  can  be  extended  to  larger  clusters  (e.g.,  the 
tetrahedron)  systematically. 

The  result  of  the  tetrahedron  treatment  in  Appendix  2 was  compared 

12 

with  that  of  the  pair  treatment  done  many  years  ago.  We  also  did  SP 
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calculations  and  compared  the  results  with  others.  These  different 
methods  compare  nicely,  and  the  relations  among  them  are  understandable. 
As  expected,  the  SP  method  can  pinpoint  the  transition  temperature 
within  the  boundary,  but  the  sum  method  cannot. 
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SECTION  b 


PHASE  BOUNDARY  AT  T - 0 

In  Appendix  5,  we  study  the  structure  of  boundaries  in  the  Cu-Au 
alloy.  The  work  In  Appendix  5 was  made  possible  by  the  previous  theore- 
tical calculation  of  the  phase  diagram  of  this  alloy  worked  out  by 
Kikuchi  and  deFontaine^’ ^ extending  the  work  by  van  Baar.^ 

During  the  work  discussed  in  Appendix  5,  it  became  evident  that 
many  properties  of  the  alloy  near  T-0  are  singular  theoretically  and 
hence  that  more  detailed  knowledge  is  needed  than  that  given  in  Refs.  15 
and  lb  about  the  behavior  near  T ■ 0 of  the  Cu^Au  phase  and  also  about 
the  phase  boundary  between  Cu^Au  and  the  disordered  phase.  Thus,  we 
did  tlie  study  near  T “ 0 and  discussed  our  results  in  Appendix  A.  The 
phase  diagram  ot  the  Cu-Au  alloy  in  Appendix  5 is  calculated  using  the 
tetrahedron  as  the  basic  cluster  of  OVM  and  also  using  the  multiatomic 
interaction  potential.  The  latter  means  that  tetraliedra  CuCuCuAu  and 
CuAuAuAu  contribute  different  four-body  interaction  potentials.  The 
four-bodv  interaction  can  induce  the  asymmetry  observed  in  experiments 
between  the  Cu-slde  and  the  Au-slde  in  the  phase  diagram. 

The  t on i - bod  v effect  is  represented  by  two  parameters  nt  and  R. 

The  range  ot  values  ot  i and  f'  in  which  the  phases  Cu^Au,  CuAu,  and 

CuAu.  .ire  stable  w.is  investigated.  The  technique  of  analysis  was  the 

17  18 

linear  pt  ograom  I ng  method,  which  had  been  used  by  Cahn  before.  * 

Ibis  part  t ’i.  wotk  wa  .lone  bv  Dr.  J.W.  Cahn  of  the  National  Bureau 
ot  Standard  , w • .'operating  in  this  project. 

In  al  nil’  a ■.  diagrams,  it  is  often  helpful  to  know  phase 

boundaries  at  V - 0.  Iliet  was  no  theory  to  treat  this  problem  before. 
Appendix  • show,  how  to  do  it.  At  a finite  temperature,  a phase 
boon  it',  is  , 1 1 ii  I a t > <1  bv  drawing  a common  tangent  to  free-oncrgv  curves 
Fj  and  F,  tor  the  two  phases  plotted  against  the  composition.  As  d i s- 
coveted  In  Appendix  A,  the  phase  boundary  between  the  disordered  phase 
anti  t hi-  ordeted  Co  (Au  phase  can  be  calculated  bv  the  following  procedure. 


15 


Configurations  of  the  tetrahedral  cluster  are  limited  to  CuCuCuCu 
and  CuCuCuAu  only  (the  rest  of  the  configurations,  e.g.,  CuCuAuAu, 
having  zero  probability  of  appearance).  Using  these  two  configurations, 
we  write  the  entropy  expressions  for  the  disordered  phase  Sj}  and  for 
the  ordered  phase  as  functions  of  the  composition.  Figure  7 of 
Appendix  4 shows  an  example.  Then  we  draw  a common  tangent  to  the 
and  curves.  The  points  of  contact  of  the  common  tangent  to  and 
give  the  two  composition  values  of  the  coexisting  phases 

The  procedure  of  finding  the  common  tangent  to  S curves  Is  equiva- 
lent to  the  following.  At  a finite  temperature,  constructing  the  com- 
mon tangent  to  the  F curves  Is  equivalent  to  finding  the  Intersection  of 

A 

two  grand  potential  G curves  for  the  two  phases  plotted  against  the 
chemical  potential  p.  When  T Is  Infinitesimally  small,  p and  G near 
the  phase  boundary  can  be  expanded  as 

p ■ p + akT  + ... 
o 

(5) 

G - G + bkT  + . . . 
o 

A 

where  p ^ and  Gq  are  common  to  the  two  coexisting  phases.  The  coef- 
ficient b Is  derived  when  the  value  of  a Is  assigned.  Then  two  b curves 
for  the  two  phases  are  plotted  against  a to  find  the  intersection,  which 
gives  the  coexisting  phases  for  the  limit  T -*■  0.  It  Is  illustrated 
in  Figure  8 of  Appendix  4. 


E 


I 
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SECTION  7 


INTERPHASE  AND  ANTIPHASE  BOUNDARIES  IN  Cu3Au 

The  f.c.c.  lattice  can  be  separated  into  four  equivalent  simple 
cubic  sublattices.  When  one  of  the  s.c.  sublattices  is  preferentially 
occupied  by  Au  atoms  and  the  other  three  s.c.  sublattices  are  equivalent, 
we  obtain  the  Cu^Au  structure  (also  called  the  Ll2  structure). 

An  APB  in  Cu^Au  is  formed  when  the  Au  atoms  on  the  left  side  of 
the  boundary  preferentially  occupy  sublattice  1 and  those  on  the  right 
side  preferentially  occupy  sublattice  2,  (See  Figures  1,  2,  and  3 of 
Appendix  5.)  When  the  left  side  of  the  boundary  is  the  Cu^Au  phase  and 
the  right  side  is  a disordered  phase,  we  call  the  boundary  the  inter- 
phase boundary  (IPB). 

APBs  and  IPBs  are  studied  in  Appendix  5 using  the  sum  method.  The 
technique  is  an  application  of  one  developed  in  Appendix  2.  We  take  a 
tetrahedron  as  the  basic  cluster  and  assign  a variable  Zn(i,j,k,jt)  to 
the  configuration  (i,j,k,ft)  of  the  tetrahedron.  The  Cu  and  Au  atoms 
are  designated  by  i = 1 and  2,  and  n indicates  the  location  of  the 
tetrahedron  relative  to  the  boundary. 

The  grand  potential 

G - F - (y1N1  + u2N2)  (6) 

for  the  entire  system  including  the  boundary  region  is  written  in  terms 
of  Zn(i,j,k,&)  and  is  minimized  (keeping  T and  y fixed)  with  respect 
to  the  Zs.  (We  can  choose  -y^  **  y = y.)  The  resulting  equations  are 
solved  for  Zn(i,j,k,£)  using  the  NIM.  The  excess  free  energy  0 attributed 
to  the  boundary  is  calculated  as  the  difference  between  G thus  calcu- 

A 

lated  and  the  value  of  G for  the  homogeneous  phase. 

The  a curves  for  the  APBs  are  given  in  Figure  6 of  Appendix  5. 

The  special  features  are  the  following: 

(1)  A 0 curve  for  a constant  y increases  monotonir.all  y as 
T decreases.  Although  in  the  phase  diagram  (Figure  5 
of  Appendix  5)  all  the  y-constant  curves  converge  to 
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the  stoichiometric  state  x ■ 1/4  as  T -*■  0,  the 
a curves  for  different  us  tend  to  different  values 
as  T *■  0.  This  is  one  of  the  singular  behaviors  of 
the  Cu.jAu  phase  near  T “ 0. 

The  lower  two  o curves  in  Figure  h of  Appendix  5 are 
for  constant  composition  x,  the  x being  the  value  for 
the  peak  (the  congruent  point)  of  the  CU3AU  phase. 

Different  from  the  p-constant  curves,  these  two 
o curves  go  through  maxima  at  about  0.6  of  the  dis- 
ordering temperature  and  reduce  to  zero  at  T - 0. 

The  general  shape  of  the  lower  two  0 curves  in 
Figure  6 of  Appendix  5 resembles  that  of  the  shear  strength 
curves  of  Nl-based  superalloya  (e.g,,  Ni'jAl).19«-0 
The  resemblance  Implies  that  the  APB  behavior  is 
an  important  factor  In  understanding  the  shear 
strength  of  these  superalloys.  The  APB  and  the  shear 
strength  are  related  through  the  creation  of  an  APB 
when  a dislocation  enters  a dispersed  ordered  preci- 
pitate within  the  disordered  matrix. 

The  disorder-CujAu  IPBs  were  calculated  and  are 
shown  in  Figures  7 and  8 of  Appendix  5.  The  note- 
worthy feature  is  that  the  iPB  is  exactly  one-half 
of  the  APB  at  the  same  point  in  the  phase  diagram. 

In  other  words,  when  the  disordered  phase  (D)  coexists 
with  the  C113AU  phase,  the  APB  inside  the  CU3AU  phase 
at  tills  point  is  made  of  two  D-CU3AU  IPBs.  This 
property,  which  is  also  seen  by  comparing  Figure  11(a) 
with  Figure  11(b)  in  Appendix  5,  means  that  the  APB 
at  this  point  is  completely  wet  with  the  disordered 
phase . 

The  density  profiles  across  the  boundary  are  shown  in 
Figures  11(a)  and  12  of  Appendix  5.  Although  the 
bulk  phase  is  C113AU  (IJ2)  structure,  there  is  a 
CuAu  (Ll0)  type  region  near  the  center  of  the  boundary. 

The  onset  of  the  Ll0  region  is  interpreted  as  a second- 
order  phase  change  from  the  behavior  of  the  plane  next 
to  the  center  as  shown  in  Figure  14  of  Appendix  S. 


SECTION  8 


GRAIN  BOUNDARY  STUDIES 

As  the  first  step  in  studying  the  structure  of  grain  boundaries, 
we  worked  on  a two-dimensional  lattice-gas  model  that  is  capable  of  pro- 
ducing gas,  liquid,  and  solid  phases  and  also  two  different  orientations 
of  the  solid  phase.  The  region  near  the  grain  boundary  is  shown  in 
Figure  1 of  Appendix  6. 

The  grain  boundary  free  energy  0 for  a constant  chemical  potential  y 

is  shown  in  Figure  3 of  Appendix  6.  Along  with  0,  we  calculated  the 

excess  entropy  S due  to  the  grain  boundary.  A remarkable  discovery  is 

that  the  S curve  is  clearly  made  of  two  portions,  as  shown  in  Figure  5 

of  Appendix  6.  For  low  temperatures,  T < 0.3  T^  (T^  being  the  melting 

temperature  for  this  y),  S is  almost  equal  to  kHn2.  For  high  temperatures, 

T > 0.6  T , S is  linear  in  -log(T  - T) , diverging  at  T . 

m mm 

This  remarkable  property  of  S,  combined  with  the  calculated  values 
of  0 for  the  grain  boundary  and  the  IPB  (the  boundary  between  the  solid 
phase  and  the  liquid  phase),  indicates  that  the  grain  boundary  is  com- 
pletely wet  at  the  melting  point. 
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Natural  iteration  method  and  boundary  free  energy* 

Ryotchi  Kikuchi 
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The  neiural  iteration  lechuique,  which  vn  proposed  by  the  author  ae  a method  of  wiving  equations  nf 
the  cluster  variation  method,  is  extended  to  the  cue  in  which  variables  are  subject  to  subsidiary  conditions 
due  to  symmetry  requirements;  the  subsidiary  conditions  are  treated  by  minor  iterations  at  each  step  of  the 
major  iteration.  The  technique  of  evaluating  the  second  order  transition  point  for  such  problems  is 
presented,  the  technique  applies,  in  general,  to  problems  for  which  the  number  of  variables  is  large.  As 
examples,  cluster  variation  methods  using  a five-point  If -shaped  duster  and  a six-point  double-square 
cluster  are  presented  for  the  two  dimensional  I sing  model.  An  improved  proof  of  the  --.l.r  product  (SP) 
expression  for  the  boundary  free  energy  is  given.  The  results  of  the  examples  are  used  in  the  SP  expression 
to  evaluate  the  boundary  free  energy  of  the  Ising  model  in  order  to  test  the  accuracy  of  the  approximations 
involved. 


I.  INTRODCUTION 

This  paper  has  dual  purposes.  One  is  to  show  in  the 
natural  iteration  (NI)  method  how  to  treat  subsidiary 
conditions  and  how  to  calculate  the  second-order  transi- 
tion point,  and  the  other  is  to  report  further  results  ob- 
tained by  the  scalar  product  expression  of  the  boundary 
free  energy.  These  two  are  closely  tied,  as  the  results 
of  the  former  are  used  in  the  calculation  of  the  latter. 

The  cluster  variation  (CV)  method1**  for  cooperative 
systems  had  been  in  use  for  many  years.  One  basic  dif- 
ficulty which  had  prevented  wide  use  of  the  method  had 
been  the  step  to  solve  simultaneous  algebraic  equations 
of  high  orders.  * This  difficulty  had  been  dissolved  by 
the  natural  iteration  method. 4 Many  applications  of  the 
new  technique  have  been  and  are  currently  worked 
out.  *-* 

The  applications  of  the  NI  method,  however,  have 
been  limited  to  either  the  pair  approximation  or  the  tet- 
rahedron approximation  (both  for  the  fee  and  the  bcc  lat- 
tices), which  do  not  need  subsidiary  conditions  on  the 
probability  variables  except  for  the  normalization  to 
unity.  As  the  application  of  the  NI  method  widens,  need 
arises  that  additional  conditions  or  the  probability  vari- 
ables are  to  be  taken  into  account  without  damaging  the 
versatility  of  the  NI  method.  It  will  be  shown  in  Secs. 

0 and  in  that  subsidiary  conditions  due  to  symmetries 
of  variables  can  be  treated  by  way  of  iterations  of  a kind 
similar  to  the  main  NI  procedure. 

One  of  the  advantages  of  the  NI  method  is  that  it  is 
unnecessary  to  carefully  choose  the  independent  vari- 
ables and  then  to  write  the  relations  among  the  depen- 
dent and  independent  variables.  When  we  work  with 
first-order  phase  transitions,  there  is  no  particular 
complication.  In  the  case  of  the  second-order  transi- 
tion, however,  the  critical  point  Te  is  to  be  calculated 
as  the  point  at  which  a certain  determinant  vanishes. 

A method  of  calculating  elements  of  the  determinant, 
without  explicitly  listing  relations  among  the  dependent 
and  independent  variables  and  hence  without  spoiling  the 
spirit  of  the  NI  technique,  is  presented  in  Sec.  II.  It 
is  believed  timely  and  useful  to  report  on  these  tech- 
niques of  handling  subsidiary  conditions  and  the  deter- 
minant for  Te  because  the  interest  in  the  CV  method  is 


growing.  *~14  The  iteration  method  recently  proposed  by 
Weeks  and  Gilmer14  for  studying  boundary  structures 
can  also  be  interpreted  as  a variation  of  the  NI  method 
with  subsidiary  conditions14 

When  the  author  wrote  a paper  on  the  scalar  product 
(SP)  expression  of  the  boundary  free  energy, 14  supple- 
menting the  work  by  Clayton  and  Woodbury,"  the  proof 
of  the  SP  expression  was  incomplete.  For  that  reason, 
the  surface  tension  o of  the  boundary  between  two  two- 
dimensional  Ising  spin  phases  was  calculated  using  the 
proposed  method  and  was  compared"  with  the  exact  re- 
sult due  to  Onsager.  Since  the  publication,  the  proof  of 
the  SP  expression  has  been  improved  and  It  is  reported 
in  Sec.  TV.  Results  of  Secs.  II  and  in  are  used  in  the 
SP  expression  to  calculate  the  boundary  free  energy  in 
Sec.  V. 

II.  THE  ^-APPROXIMATION 
A.  Free  energy  and  its  minimization 

We  present  discussions  on  the  NI  method  in  this  sec- 
tion using  the  W-approximation  of  the  two-dimensional 
Ising  net  as  an  example.  This  approximation  uses  a 
five-point  cluster,  A-B-C-D-E  of  Table  IA,  as  the 
basic  cluster.  The  degeneracy  factor  ft,  for  this  case 
was  shown  in  Table  I of  Ref.  2 and  is  reproduced  in  Ta- 
ble IB.  Reference  2 is  to  be  consulted  for  the  meaning 
of  the  parentheses  notation  in  R*.  Expression  Rr  is 
used  in  writing  the  entropy  of  the  system  as 

S=  felnfty  . (2. 1) 

The  plus  and  minus  spins  are  designated  by  i - 1 and  2, 
respectively.  The  probability  of  finding  a spin  configu- 
ration i-j-k-l-m  on  A-B-C-D-E  points  of  the  cluster 
is  written  as  ulmm,  as  shown  in  Table  IC.  Other  vari- 
ables, v,  z and  y,  are  also  defined  in  Table  IC.  Using 
these  variables,  the  entropy  (2. 1)  for  a system  of  A* 
lattice  points  can  be  written  explicitly  as 

+ ,v„)  , (2.2) 

where  the  £ operator  is  defined  as 


1 


f*i 
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TABLE  I.  The  *■  c Uieler. 
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Each  summation  In  (S.  2)  is  to  be  done  over  the  values  1 
and  2 of  the  subscripts  In  the  summand.  The  two  v 
terms  In  (2.  2)  are  equal,  but  are  written  separately  In 
order  to  make  the  subsequent  NI  formulation  symmetric. 

In  a system  of  S lattice  points,  the  total  number  of  H’- 
clusters  we  are  Interested  In,  i.e.,  those  lying  parallel 
to  each  other,  la  2V.  The  energy  of  the  system  Is  then 
written  as 

E -2A- , (j.4) 

where  la  the  energy  per  H’-ctuster  and  is  written, 

In  turn,  in  terms  of  pair-wise  energy  tlt  as 

•wie'Htii+ta+tu+ti,)  . (2.5) 

We  define  «(<  as 

when  i->,  (J  e) 

<<!*♦  « , when  i*j  . 

Our  program  is  to  minimlie  the  free  energy  with  re- 
spect to  » UUa’i.  The  a'  variables  obey  not  only  the 
normalisation: 

E "i/*i«  * 1 . (2. 7) 

i,i,*, i, « 


IlM  “ I]  mi tUm  " E 


(2.8) 


The  mirror  symmetry  is  easy  to  take  care  of  because 
the  Nl  treatment,  to  be  described  in  detail  below,  keeps 
this  symmetry  intact  through  the  iterations  tf  the  initial 
input  it’s  satisfy  the  symmetry. 

The  translational  symmetry  is  taken  into  account  by 
adding  the  following  Lagrange  terms  to  the  free  energy 
expression: 

^ • * XI  °ISI  X ^“’imi.  * •*'«li»|)  » 

m 

* • (2.10) 

t, 

We  can  now  prove  that  the  Lagrange  multipliers  a's 
obey  the  following  symmetry  relations: 

«4JM  ■ * *1*1  • (Ml) 

The  proof  is  the  following.  We  form 
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(2.12) 

Ws  see  thst  the  first  summand  in  vanishes  when  the 
mirror  symmetry  (2.8)  is  satisfied)  this  allows  us  to 
choose  + alkJ,  ■ 0,  which  is  (2. 11).  This  relation 
(2.  11)  was  also  proved  in  (U.9)  of  Ref.  2 in  general 
notation. 

Using  relations  obtained  so  far  in  this  section,  we 
can  write  the  free  energy  K - TS  as 

jiF  fiE  _S_ 

N S~  kN 


-2  E (<* 

Ml* 


Ml  < 


where 


(2. 13) 


and  we  use  E of  (2. 4)  and  S/kS  of  (2.  2)  in  (2. 13).  We 
then  minimise  F in  (2. 13)  with  respect  to  In 

order  to  keep  the  symmetry,  we  use 

«'o»i  * E *tiMm  > 

m 

and 


«>«.*  E*'«i»i.  . (2.14) 

i 

for  e,  „ and  tn  (2.  2).  The  differentiation  of  (2. 131 
leads  to 


but  also  two  symmetry  requirements.  One  of  them  is 
the  "mirror”  symmetry: 


where 


t2.  15a) 


tlrtlK  ' ailUl  » 

and  the  oil.  r is  the  "translational”  symmetry: 


(2. 8)  ♦ I ln(t  0»,  t MtK *UK/V<I)  . (2. 15b) 

We  have  Introduced  reJ5»„.for  the  sake  of  convenience  in 
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the  following  discussions. 

Equation  (2. 15)  is  the  "superposition"  relation*  which 
is  ready  to  be  used  in  the  NI  treatment.  Except  for  a’ a, 
which  will  be  discussed  in  the  next  subsection,  the  Nl 
starts  with  an  input  set  {ivl>t,a },  and  calculates  the 
right-hand  side  of  (2. 15a)  to  obtain  the  output  set 
{“'imi-}*  The  input  t/’s  are  derived  from  w’a  using 
(2.  14),  and  the  input  x’s  and  y's  are  from 

zltm~  *’(/*!«  I 

hi 

„ (*• 1«> 

y»i  - 2-.*  i u • 

* 

In  order  to  start  the  NI  procedure,  we  assign  the  val- 
ue of  kT/i  and  initial  guess  of  Instead  of  speci- 

fying all  32  ic’s  independently,  it  is  sufficient  to  start 
with  giving  probabilities  of  a single  spin  x,  and  x,  (x, 

+ x,  = 1)  and  then  write  w’a  as  a Bragg- Williams-type 
product: 

(2.17) 

This  initial  assignment  works  even  for  low  tempera- 
tures. 

We  call  the  iteration  step  of  going  from  the  input  set 

to  the  output  set  which  is  the  next  input 

set,  the  “major"  iteration  step.  Convergence  of  the 
major  iteration  is  tested  by  the  following  sum  for  each 
step: 

= X. 1 )n,*’i4»i«  ~ *n“’u»i«il  • (2.18) 

A typical  example  of  AmJ  as  a function  of  the  major  It- 
eration step  is  shown  in  Fig.  1.  It  is  better  to  use  a 
logarithm  In  defining  4a,  as  in  (2. 18)  than  using  a sim- 
ple sum  £1  futim-  J u8ed  1”  Ref.  4,  because  some 
of  the  ic's  become  very  small  for  low  temperatures  and 
do  not  contribute  to  I w-  u>| , although  they  do  to  (2. 18). 
We  need  accurate  digits  of  these  small  tv’s,  particularly 
in  calculating  the  boundary  free  energy  using  the  scalar 
product  expression  in  Sec.  V. 

In  the  present  problem,  Fig.  1 shows  that  4^  be- 
haves in  the  same  way  as  was  reported  in  Ref.  4 and 
logh^i  decreases  linearly  as  the  iteration  progresses. 
It  is  safe  to  use  4al  = 10*1  as  the  criterion  of  the  con- 
vergence. 

We  can  prove  that  when  the  Iteration  has  converged, 
the  expression  (2. 13)  reduces  to 

» = F/N  . (2. 10) 

B.  The  minor  iterations 

The  remaining  question  In  the  present  case  is  how  to 
handle  a’s  which  are  to  be  determined  so  that  the 
translational  symmetry  (2.0)  is  satisfied.  When  we 
substitute  (2. 15a)  in  (2.9),  we  obtain 

ai4»i  = i Ms*  (2. 20a) 

who  i e 

exp(o  ■itjft)  » (2*  20b) 


as 

Our  procedure  of  satisfying  (2. 20)  is  to  do  a series  of 
“minor”  iterations  to  determine  a’s  at  each  "major" 
iteration  step. 

In  the  method  we  propose,  we  use  input  a’s  on  the 
right-hand  side  of  (2.20a)  and  evaluate  aIJtl  on  the  left- 
hand  side  as  the  output  of  a minor  iteration  step.  This 
output  is  used  as  the  next  input  on  the  right-hand  side. 
This  technique  for  the  minor  iteration  is  of  the  same 
spirit  as  the  NI  major  iteration  steps  and,  although  we 
have  not  analytically  proved  convergence  of  the  minor 
iterations,  they  do  converge  in  all  similar  cases  we 
have  worked  so  far.  Speaking  in  qualitative  terms,  the 
output  a,,,,  in  (2.  20)  is  derived  as  a weighted  average 
of  input  a’s  on  the  right-hand  side;  the  averaging  pro- 
cess is  Interpreted  as  helping  convergence. 

At  each  major  iteration  step,  the  minor  iterations 
are  tested  using 

Aw=  21  I aij»i  ~ ®<4»i  I • (2.21) 

14*1 

When  becomes  less  than  a criterion  value  Aa<r>#, 
the  Iteration  is  judged  converged.  The  number  of  mi- 
nor iterations  for  one  major  iteration  step  depends  on 
AuraC  and  also  on  the  progress  of  the  major  iterations, 
and  gradually  reduces  to  one  as  the  major  iteration  ap- 
proaches its  convergence.  A typical  example  of  the  be- 
havior of  the  minor  iteration  is  shown  in  Fig.  2. 

As  an  alternative  method,  the  simultaneous  nonlinear 
equations  In  (2. 20)  for  a’s  can  be  solved  by  the  Newton- 


MAJOR  ITERATION  STEP 

FIG.  1.  Convergence  pattern  of  the  major  Iterations. 
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Rapnson  (N-R)  method.  The  convergence  is  expected 
to  be  good  because  both  positive  and  negative  values  of 
a’s  are  allowed  and  there  is  no  danger  of  overshooting 
into  un physical  regions.  Different  from  the  method  we 
are  proposing  in  this  section,  however,  the  N-R  pro- 
cedure needs  evaluation  of  derivatives  and  an  inversion 
of  a matrix.  In  the  present  example,  the  N-R  method 
is  not  so  bad  because  the  number  of  independent  a’s  is 
only  six  m the  ordered  phase.  For  problems  in  which 
the  number  of  independent  a’s  is  larger,  the  method  we 
are  proposing  is  faster. 


C.  Determination  of  Tc 

For  the  sake  of  convenience,  we  call  a plus  spin  and 
a minus  spin  complements,  and  write  the  complement 
of  the  subscript  i as  i'  (f  ♦#'  * 3).  The  variables 
and  My become  identical  in  the  disordered  phase 
and  we  define  their  difference  as  the  long-range  order 
parameter  4*»  (For  ( we  use  a different  subscript  sys- 
tem.) Because  of  the  mirror  symmetry,  the  number  of 
long-range  order  parameters  is  ten.  An  example  is 

* "'mu  ~ "11*11  • (2. 22a) 

The  list  of  subscripts  of  the  first  term  m’s  is  (11111), 
(11112),  (11121),  (11211),  (11122),  (11212),  (12112), 
(21112),  (11221),  and  (12121)  in  the  order  of  p = 1, 

2, . . . , 10.  It  is  convenient  to  rewrite  (2.  22a)  as 

&im  h~  ("'mu  ~ Rfeitti)  * 0.  (2.22b) 

Beside  the  differences  of  m’s  differences  of  a’s  also 
vanish  in  the  disordered  phase  and  can  be  defined  as  the 
long-range  order  parameters  4„(p  = 11,. .. , 14).  There 
are  four  of  them: 


^11  * (|1  — (<*uu  — <T«**i)  = 0 
tfll  5 s|l  - (<*Hi|  — <*UI|)  - 0 

R\i  * tu  ~ (anu  - ®mii)  * 0 
R 14  * i|«  — (®l*|I  “ - 0 


(2.  23) 


The  second-order  transition  point  T,  is  determined 
as  the  |X)int  at  which  the  following  14  x 14  determinant 
vanishes: 


Det 


= 0. 


(2.  24) 


Since  we  are  working  at  T^,  the  long-range  order  pa- 
rameters 4„  are  small.  Thus  we  also  write,  along  with 
(’s  interchangeably, 


Am 

Aar 


Dtl  i 


,*  ", 


■ 


I/Mi 

11*1  * aUkl  ~ ®l>/<»'|< 


(2.25) 


The  expression  gm  as  defined  in  (2.  22b)  and  (2.  23)  can 
be  expanded  in  terms  of  these  A’s.  in  other  words,  we 
can  * rite  in  principle 

it 


’ I'm  (*  • 


(2.26) 


The  coefficient  <■„,  is  8#„/84,„  which  we  want  in  (2.24). 
In  order  to  evaluate  r„„  we  do  not  need  to  derive  the 
expression  (2.26)  analytically  but  can  use  the  following 
trick.  We  simply  choose  1 and  4,  = 0 for  « and 
using  substitutions  (which  we  explain  below)  of  many 
steps,  we  finally  calculate#,,.  However  complicated 
the  substitutional  steps  may  be,  the  value  of  thus  cal- 
culated is  equal  to  c„. 

We  show  the  substitutional  steps  using  of  (2.22b) 
as  an  example.  For  the  quantity  An  lilt,  in  the  paren- 
these  in (2. 22b),  weusetherelationderlvedfrom(2. 15) 


Air, 


/Aw 


V0> 


+ Aar,,!,  * A a, 


.)  • (2. 


27a) 


In  deriving  this  expression  from  (2.  15a),  we  left  out 
A(>0)  because  \f)  changes  smoothly  across  Tc  and  hence 
I 8U/3)/87-|  « I »(,JdT\  at  T The  first  term  in  the 
parentheses  in  (2.27a)  is  calculated  from  (2. 15b)  as, 
again  leaving  out  A0, 

Aim*01 
111 


(2.  27b) 


-“jiiui  = I/^luu  + + A*nt  _ A.vu  \ 

"’uiu  2\t>u„  ty,w  *ui  v„  /’ 

where  we  use  differentiations  of  (2.  14)  and  (2. 16)  for 
Ar’s,  Ar’s,  and  Ay’s  to  write  them  as  linear  combina- 
tions of  Ax’s.  For  Aatlu  and  Aalm  in  (2.27a),  we  use 
the  differentiation  of  (2.  20) 


Aar,,,,  = 

,<w  »\  «*..,*,  s, 


i Jkl 


)• 


(2.  28a) 


^S/r.i/* i = H (a»-i/*.  + R'-ij» iA«-ij*)e*P(a«iM)  , 


Aa,„() exp(a-iw)  . 

(2.  28b) 

The  transition  point  thus  determined  from  the  vanish- 
ing of  the  determinant  in  (2.  24)  is 

*T«A  = 2.  36483  , (2.29) 

which  is  the  same  as  the  value  reported  in  Ref.  2 ex- 
cept for  the  last  digit  3 which  was  reported  as  0 previ- 
ously. 

III.  THE  DOUBLE  SQUARE  APPROXIMATION 
A.  Free  energy  and  its  minimization 

As  was  discussed  at  length  in  Ref.  2,  the  tl-approxi- 
mation  treatment  in  Sec.  II  looks  at  the  square  net  from 
the  (11]  direction,  and  is  appropriate  in  calculating  the 
(11)  boundary  free  energy.  When  we  treat  the  lattice 
from  the  [10]  direction,  the  basic  cluster  which  is  one 
step  larger  than  the  square  is  a double  square  as  shown 
by  A- B-C- D—E—F  in  Table  I LA.  This  cluster  is  ap- 
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TABLE  II.  The  double-sryiare  cluster. 


propriate  in  calculating  the  (10>  boundary.  The  ilDS  fac- 
tor for  this  case  is  also  shown  in  Table  UB.  This  clus- 
ter was  studied  independently  by  Allegra  and  Delise  in 
Politechnico  di  Milano. 19 

The  probability  variables  are  defined  in  Table  llC. 

As  in  Sec.  II,  i = 1 and  2 denote  plus  and  minus  spins, 
respectively.  Using  these  variables  and  the  (1DS  factor, 
we  can  write  the  entropy  S for  a system  of  N lattice 
points  as 

,)  + o + i 

- £(  y J + y*J + £<e(  y>  i>  + yfJ)  • 

(3.1) 

where  the  £ operator  is  defined  in  (2. 3)  and  the  sum- 
mation of  each  term  goes  over  the  values  1 and  2 of  the 
subscripts  in  the  summand. 

In  a system  of  N lattice  points,  the  number  of  double- 
square  clusters  is  N.  Thus  the  energy  of  the  system 
is  written  as 

* =^X/er/»r«M  mm  t (3.2a) 

where  « is  now 

i(«<i*-C»i  + <*,,)♦  i (<«+ t/j +€*»,+ CjJ  • (3.2b) 

The  pair  energy  c<y  is  defined  in  (2.  6). 

We  now  examine  subsidiary  conditions  on  w’a.  The 
normalization  is 

Z • (3.3) 


There  are  two  “mirror”  symmetries: 

and 

- Wftlkmm  • (3.4) 

These  two  symmetries  in  (3.4)  can  be  kept  through  the 
iterations,  if  they  are  satisfied  in  the  starting  input  «•’ s. 
There  are  two  more  symmetries.  One  is  the  “transla- 
tional” symmetry: 

= X (3.  5) 

at,  n *»,  * 

and  the  other  is  the  “rotational”  symmetry  of  f’s; 

vn»i  = vttn  • (3.  6) 

The  last  one  is  to  guarantee  that  the  nrobability  for  the 
spin  configuration 


is  equal  to  that  of 


and  hence  the  isotropy  of  the  system. 

The  translational  symmetry  is  taken  into  account  by 
Lagrange  terms  similar  to  £„  in  (2. 10): 

£«  = S (®(i*l  + aau*|)  u’ttklmn  l (3.7) 

where  we  have  used  the  symmetry  of  a’ s: 

atjt,  = -atllJ . (3.8a) 

This  relation  can  be  proved  by  an  argument  similar  to 
(2. 12).  We  may  note  that  a’s  obey  another  symmetry 
relation: 

aiMi  ~ ann  » (3. 8b) 

which  corresponds  to  the  second  symmetry  of  w’s  in 
(3. 4).  Because  of  these  two  symmetry  requirements, 
the  Independent  a’s  are  three:  auu,  ann,  and  altu. 
Six  a’s  vanish. 

The  rotation  symmetry  (3.6)  needs  the  following 
Lagrange  terms: 

= X (t* tiki + Ymnkt)  wi/klmn  f (3.  9) 

U»l 

where  y’s  satisfy  the  symmetry  relations: 

Yimi  ~ ~ You  I > 
and 

Yi/h  =Yjut  . (3.10) 

Because  of  these  relations,  12  y’s  vanish  and  there  is 
only  one  independent 

Yutt~  Ytui  - —Yiiu  - — V-2U  • (3.11) 

The  free  energy  expression  corresponding  to  (2.  13) 
is 
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QF  or  J 

~n  * ~n  ~ On  ~ + + i + 

t 

+ X/}  ^1  — ^ i (3.  12) 

where  we  use  E of  (3.  2)  and  S/kN  of  (3. 1).  When  we 
minimize  (3. 12)  with  respect  to  we  obtain 

Inx’uuw.  = *0  + lme}*i,„+  ®u»i  + a«*i + r<i*i + V-wi  . 

(3. 13a) 

where 

* i ln(«q,»i  ~ i ln(  IWji  J'J  • 

(3. 13b) 

Except  for  the  minor  iterations  to  be  shown  in  the 
next  subsection,  the  natural  iteration  step  starts  with  an 
input  set  {wt„lmn\,  calculates  the  output  set  {u>IJilm„} 
using  (3. 13),  and  then  uses  w’s  as  the  next  Input.  The 
input  r’s  are  calculated  from  w’s  in  (3.  5),  and  the  input 
z’s  and  y’s  are  from  the  following 

*l»m~  ^2  " ijtlmn  > 

J.l.n 

(3.14) 

yin  ~ • 

m 

We  will  call  the  iteration  step  from  the  input  set 
{“'lini.ii}  to  the  next  input  set  the  major  iteration.  We 
can  use  an  expression  similar  to  (2. 18)  as  the  test  value 
naturally  with  the  six  subscripts  ( ijklmn ) in  the 
present  section  rather  than  five  in  (2. 18).  The  con- 
vergence behavior  in  the  present  section  is  similar  to 
Fig.  1,  but  slightly  faster. 

When  the  iteration  converges,  the  value  of  X is  again 
equal  to  F/N  as  was  the  case  in  (2. 19). 

B.  The  minor  iterations 

For  each  major  iteration  step,  the  output  iti’s  have  to 
satisfy  the  "translational”  and  the  “rotational”  sym- 
metry relations  (3.  5)  and  (3. 6).  We  have  to  solve  the 
following  two  sets  of  equations  for  three  a’s  and  yu„: 

aun  = i l|,('Sj¥1n*i/Sptij»|)  > (3. 15a) 


Stf,  ll»l  = wnnl)m  exP(a«mli  + V«ii,l>)  » 

MR 

Sd.UH  = 5Z  “{i*! mu exp( + y*.,,)  , 


(3. 15b) 


(3. 16)  to  evaluate  the  output  set  {a,  y}  and  repeat  tin- 
cycle  until  the  iteration  converges.  At  each  major  it- 
eration step,  the  minor  iterations  are  tested  using 
similar  to  (2.  21)  but  including  the  y term.  The  minor 
iterations  behave  somewhat  similar  to  Fig.  2,  but  re- 
duce faster  to  one;  in  the  example  of  the  same  temper- 
ature, the  points  reduce  to  one  in  Fig.  2 around  80-100 
major  iterations. 

C.  Determination  of  Tc 

We  can  determine  Tc  using  the  method  similar  as 
Sec.  1IC.  Corresponding  to  (2.  22),  there  are  ten  equa- 
tions of  the  form 

Ku*  — (“'l/HIMB  ~ “Vl'H'I'.'n')  - 0 . (3.  17) 

The  equations  corresponding  to  (2. 23)  reduce  to  two: 

A'n  s in  - (<*111*  “ ami)  » 

^uMi* -(®ii**  - <***u)  > ^ ^ 

and  the  determinant  (2. 24)  now  becomes  12  x 12.  The 
"rotational”  symmetry  parameter  y„12  does  not  con- 
tribute to  this  determinant  because  ymt  - yml  = 0 even 
in  the  ordered  phase  as  is  seen  in  (3. 11).  The  value  of 
Tc  we  obtain  from  vanishing  of  the  determinant  is 

kTc/(  = 2. 37619  . (3.19) 

It  is  of  interest  to  see  the  calculation  for  which  the 
“rotation”  symmetry  is  not  imposed.  We  find  the  dif- 
ference is  relatively  small.  The  value  of  kTje  for  this 
case  is  almost  the  same  as  (3. 19)  except  that  the  last 
digit  9 is  replaced  by  3.  As  far  as  the  values  of  v's  are 
concerned,  we  have  the  following  example.  At  kT/t 
= 2.  38  which  is  just  above  Tc,  the  disordered  phase  val- 
ues are 

iijin  = — 0.  031227  , (3.  20a) 

when  the  rotation  symmetry  is  imposed,  and 

fji**  = 031373  , 

ti,»u=  0.031075  , (3.20b) 

when  it  is  not. 

When  we  compare  Tc  in  (3. 19)  for  the  double-square 
case  with  (2. 29)  for  the  W approximation,  we  find  the 
W- value  is  about  0.  5%  lower.  Based  on  the  generally 
accepted  rule,  this  means  that  the  W -approximation  is 
slightly  better  than  the  D-S  approximation.  This  dif- 
ference is  information  which  was  not  available  when 
Table  I of  Ref.  2 was  written. 


Vim  = j ln(/iy/lljj)  , (3.16a) 


Rn  ~ ^2  "'«m**ll  eXP(®UU  + 0 nnZ2  + y«mt*)  » 

Rd~  ^2  M'i*i*«nexP(ai*i*  * ° mi*  + Vimi*)  • (3*  !6b) 

«n 

In  the  iteration  technique  we  are  proposing,  the  input 
set  {a,y\  is  used  on  the  right-hand  side  of  (3. 15)  and 


D.  Comment  on  the  double-square  approximation 

The  flDS  expression  in  Table  II  is  the  case  C2  of  Table 
I in  Ref.  2,  and  treats  the  two-dimensional  net  layer 
wise.  Explaining  in  detail,  this  flDS  is  the  number  of 
ways  the  lattice  can  be  constructed  in  such  a way  that 
all  the  horizontally  placed  double  squares  (like  A-B- 
C-D-E-F  in  Table  III)  have  the  assigned  distribution 
of  spin  configurations.  It  is  important  to  note  that  this 
f!os  does  not  pay  attention  to  vertical  v placed  double 
squares  (like  A-B-G-H-D-C  in  Tab  e III). 
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When  the  interaction  energy  la  limited  to  the  nearest 
neighbor  A-K  and  the  second  neighbor  A-l>,  then  this 
l)u,  expression  is  valid.  However,  when  the  third 
neighbor  (like  A-K)  UUeractlon  is  to  be  taken  into  ac- 
count, the  i>0l  expression  in  Table  11  (alia.  The  rea- 
son la  that  the  distribution  of  the  vertical  third  neighbor 
A-C-  ia  not  taken  into  account  in  the  theory.  In  Table 
V of  his  article,’  Hurley  applied  the  Up,  exprescion  in 
Table  11  to  the  problem  in  which  the  third  neighbor  in- 
teraction cornea  into  play.  Hia  example,  therelore,  ia 
an  example  in  which  the  method  was  nut  expected  to 
work,  aa  Burley  agreed  with  the  author  privately  later. 

11  the  third  neighbor  (like  A-C)  and  the  fourth  neigh- 
bor (like  A-H)  interactions  are  to  be  included,  the  A, 
expression  of  the  U factor  in  Table  I of  Ref.  2 should 
be  used.  It  ia  reproduced  aa  fila*(A|)  in  Table  111.  This 
U treats  both  hurl  emit  a I and  vertical  double  squares 
equally,  and  ia  applicable  to  the  third  and  fourth  neigh- 
bor interactions. 

Aa  was  listed  ui  Table  I of  Kef.  2,  however,  if  the  in- 
teraction is  nearest  neighbor  only,  then  Up,  of  Table  11 
gives  a better  result  than  Utat04,)of  Table  III.  This 
indicates  that  although  the  latter  works  with  the  third 
and  fourth  neighbor  interactions,  the  approximation  may 
not  be  as  good  aa  expected.  We  certainly  expect  a good 
approximation  when  a 3x3  cluster  made  of  circled 
points  in  Table  UI  is  used.  For  thla  case  the  U factor 
Is  given  in  A,  of  Table  1 in  Ref.  2 and  is  reproduced  in 
Table  III  here  also. 

I lie  same  kind  of  argument  about  the  horizontal  and 
vertical  holds  for  the  If  approximation  of  Sec.  11. 


IV  PROOF  OF  THE  BOUNDARY  FREE  ENERGY 
EXPRESSION 

l>.  Ref.  IB  the  scalar  product  expression  for  the 

I mm  dary  excess  free  energy  was  derived.  The  proof 
of  l ie  expression,  however,  contained  a weakness.  In 
Ihf  section  we  present  an  improvement  of  the  proof. 

II  v.  ill  be  done  using  a three-dimensional  nomenclature. 
I ll  I,  we  briefly  summarize  Sec.  II  of  Ref.  16. 


We  consider  an  »th  crystalline  plane  parallel  to  llie 
boundary  and  designate  by  i>,  a configuration  ol  the  en- 
tire plane.  The  probability  that  the  ith  plane  takes  the 
configuration  i>,  is  written  as  />(i',).  We  define  the 
quantities'. 

A'i(f|)  “ I /’(*'*)  |‘ 71  exp|  i »<*»<•  ,)1  , 

which  can  be  regarded  aa  a column  vector  g,  and  a row 
vector  hlt  whose  components  are  indexed  by  e(.  The 
uKe,)  takes  care  of  the  continuity  condition  and  vanishes 
when  the  system  Is  homogeneous.  Because  of  the  nor- 
malization of  the  probability  distribution  />(r(),  the 
scalar  product  of  g,  and  h,  is  unity: 

h,  * g,  • 1 . (4.2) 

in  Ref.  IB  we  showed  that  g,  and  h,  obey  the  recur- 
rence relations: 

•i  ■ rxp(dX(,i/i)  P • gi.t  , (4.3a) 

h,  -exp(d*M/i)h,.|*  P , (4.3b) 

where  P is  the  transfer  matrix  and  depends  on  the  in- 
teraction energies.  The  parameter  X,.,/,  Is  related  to 
the  free  energy  F of  the  entire  system  as 

A"  * ^ X(.(/1  . (4.  4) 

When  the  system  is  homogeneous,  g,  * Si.i  «,0> 
that  (4. 3a)  reduces  to 

*“*-*xp(dx‘#,)P- f“"  . (4.5) 

This  Is  the  eigenvalue  equation.  The  excess  free  ener- 
gy Aa  attributed  to  the  boundary  (A  is  the  cross-sec- 
tional area  of  the  system)  is  written  as 

1 

Ao- llm  ('i,i/i - X10’)  , (4.tl) 

when  the  boundary  region  lies  tie  tween  i * - hi  and  t hi. 

From  now  on,  we  deviate  from  Ref.  16.  In  order  to 
calculate  a,  we  consider  a region  between  i « + hi  and 
♦ 3m  outside  of  the  boundary,  in  order  to  compare  with 
the  boundary  region.  We  ussume  w large  enough  so 
that  the  region  3ni  Is  sufficiently  close  to  the 

homogeneous  phase.  Thla  assumption  allows  ustowrite 

1.  (4-7) 

which  is  a modification  of  (4.2).  We  can  also  write 

g„  « expOt2m\'01)  P*"  • g,.  , (4.8) 

which  is  obtained  by  applying  (4.  3a)  or  (4.  5)  2m  times. 

The  boundary  excess  free  energy  can  be  derived  as 
follows.  Operating  (4.  3b)  between  - m and  w,  we  derive 

h.,.-oxp(-0  x,.t/1)h,  •?•*".  (4.9) 

' i •*  net  ' 

We  form  the  scalar  product  ol  4.9)  and  (4.8',  and  use 
(4.  6)  and  (4.  7)  to  arrive  at 

limit.,,*  g„  * exp(-  flAa)  . (4.  10) 

doing  buck  to  />(r4)  in  (4.  1),  \m  can  write  (4.  10)  us 
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exiH-  » (4.11) 

¥ 

where  the  subscript*  I and  U denote  the  first  and  second 
phases  which  meet  at  the  boundary. 

Although  this  proof  is  not  written  In  mathematical 
rigor,  and  still  contains  some  implicit,  albeit  reason- 
able, assumptions  how  h„  and  f . behave  away  trom  the 
boundary,  it  is  an  improvement  over  the  proof  in  Ref. 
16,  and  is  applicable  to  two-dimensional  boundaries 
for  which  the  location  of  the  boundary  cannot  be  speci- 
fied,10  as  well  as  to  three-dimensional  boundaries. 

V.  CALCULATIONS  OF  THE  BOUNDARY  FREE 
ENERGY 

When  we  apply  the  scalar  product  expression  (4.  11) 
to  the  present  problem,  we  can  write  for  the  (10) 
boundary 


exp(-«<T,1O)(/0)-  £ « {*.,*}(/>!  {*l/*}/>ll{*iy*})  , 

(5.1) 

where  a is  the  lattice  constant  and  n is  the  number  of 
lattice  points  in  a line  parallel  to  the  boundary.  The 
quantity  on  the  right-hand  side  is  the  probability 
for  a three-point  cluster  as  is  defined  in  Table  II. 

The  sum  over  all  configurations  v„  in  (4. 11)  is  now 
written  as  the  sum  over  the  set  {*,,,}  with  the  weight 
factor  This  factor  is  the  number  of  different 

ways  the  line  (composed  of  n lattice  points)  can  take 
the  configuration  specified  by  the  distribution  and 

is  written  following  the  cluster-variation  method'  for  a 
one-dimensional  system  as 


o{*i,*}= no.*,)  • / n («*i/»)  i . 

t.i  ' t.j.i 


(5.2) 


In  accordance  with  the  sum  over  in  (5. 1),  the 
probability  factor  />,(  e)  in  (4. 11)  is  changed  into 
/’t  {*u*}»  which  is  the  probability  that  a line  of  ti  lattice 
points  Imbedded  in  the  equilibrium  homogeneous  bulk 
phase  I takes  a configuration  specified  by  the  set 
The  conditional  probability  that  a spin  k Is  found  next  to 
a pair  of  spins  i-j  is  in  the  phase  I.  Making 

use  of  this  property,  we  can  write  />,  as 


II  (*i5»)**(wiii»)/ri(.V(5,)**(n.v0)  , (5.3) 

i.).>  I *.) 

where  a double  asterisk  is  the  I-ortran  notation  meaning 
"raised  to  the  power  of.  ” In  this  expression,  zjJJ  with 
a superscript  (I)  is  the  value  In  the  equilibrium  homo- 
geneous bulk  phase,  and  is  to  be  distinguished  from 
which  is  the  summing  variable  In  (5. 1).  The  other  fac- 
tor />,,  in  (5. 1)  can  be  written  by  changing  (I)  in  (5.  3) 
into  (II). 


We  are  now  ready  to  calculate  a from  (5. 1).  Using 
the  standard  technique  of  statistical  mechanics,  we  re- 
place the  summation  in  (5. 1)  by  the  maximum  term. 
For  that  purpose  we  write  from  (5. 1) 

min[£  £(*(,»)-  £l‘(y,y) 

t.J 


where  £ is  the  operator  defined  in  (2.3)  and 


(5.4) 


(5,< 

The  A/I  terms  in  (5. 4)  come  from  the  normalization  of 
when  (5.4)  is  minimized  it  reduces  to 


oiwa0=Xfi.  (5.6) 

We  note  that  the  right-hand  side  of  (5.4)  is  of  the 
form  of  a free  energy  of  a one-dimensional  system  in 
which  lny,,  and  lni(>»  play  the  role  of  energy  (times  d). 
Since  a one-dimensional  system  has  only  the  disordered 
phase,  z’s  and  v’s  obey  the  symmetry  *nt  = *ui>  etc. 
Minimization  of  (5.4)  leads  to  an  eigenvalue  equation, 
e'**  being  the  eigenvalue.  It  is  solved  as 


e 


-a* 


ijim.iiu 

2b„ 


(5.  7) 


As  a check  of  this  equation,  we  examine  the  high- 
temperature  case  when  the  homogeneous  bulk  phase  is 
disordered.  In  such  a case,  y(y  and  i tJt  of  (5.  5)  reduce 
to 

* vi}» *>  yJJ"  , and  2lit  = zJJi  = * j“>  , (5.  8) 

and  we  can  show  that  e“a  of  (5.  7)  reduces  to  unity  and 
hence  o(10)  = 0. 

The  boundary  excess  free  energy  o(I0)  for  the  (10) 
boundary  is  then  calculated  from  (5.5),  (5.6),  and  (5.7) 
using  the  results  of  the  equilibrium  homogeneous  bulk 
phase  of  Sec.  HI.  The  dimensionless  quantity  o(l0)a/« 
is  plotted  In  Fig.  3 by  a solid  curve  marked  as  double 
square.  It  is  compared  with  the  exact  result  of 
Onsager*'  and  the  results  based  on  the  pair  and  the 
square  approximations  which  are  reported  previously." 

The  boundary  free  energy  o(ll)  of  the  (11)  boundary 
can  be  calculated  using  (5.  5)  and  (5.  7)  again,  but  re- 
placing (5. 6)  by 

°<ii)  ■Ti  *0  , (5.9) 

since  the  distance  between  two  neighboring  lattice  points 
on  a (11)  line  is  vT  a.  For  this  boundary,  y(>  and  tln 
In  (5.5)  are  those  defined  in  Table  I.  Section  □ is  used 
to  calculate  o<n)  through  (5.  5),  (5.7),  and  (5. 9).  The 
result  is  shown  by  the  broken  curve  marked  with  IV  in 
Fig.  3|  It  is  to  be  compared  with  the  exact  result  by 
Fisher  and  Ferdinand**  and  also  with  the  angle  approxi- 
mation result"  which  are  also  plotted  in  Fig.  3 with 
broken  curves. 


The  curves  in  Fig.  3 support  the  correctness  of  the 
scalar  product  expression  of  the  boundary  free  energy 
as  well  as  effectiveness  of  the  IV-  and  double-square- 
approximations  of  Secs.  II  and  HI. 


VI.  SUMMARY  AND  CONCLUSION 

This  paper  can  be  regarded  as  the  follow  on  of  or  a 
supplement  to  the  several  papers  published  previously, 


J.  Chem  Phy».,  Vol.  65.  No.  It,  1 December  1976 


Ryoichi  Kikuchi:  Natural  iteration  method  and  boundary  free  energy 


4663 


I 


klU 

\ 1C.  I • udary  excess  free  energies  rr<10>  and  tr<l (>  for  the 
two  d'nien.  : >nal  Islng  model.  The  - 'lid  curves  are  for  o<ln> 
and  tl'0  broken  curves  are  o<M>.  1 >.e  kinds  of  approximation 
used  *n  evaluation  are  indicated  in  ihe  figure. 

the  paper*  on  the  improvement  of  the  cluster  variation 
method,  the  natural  iteration  paper,4  and  the  boundary 
free  energy  papers,  particularly  the  latter  two,  and 
answers  questions  left  unanswered  or  having  come  up 
after  the  publications,  in  order  to  pave  a way  for  fur- 
ther applications  of  the  teehnioues. 

In  computing  the  equilibrium  state  of  a cooperative 
system  using  the  cluster  variation  method,*  it  is  nec- 
essary in  general  except  for  simple  cases  that  subsid- 
iary conditions  are  to  be  satisfied  among  variables. 
Sections  □ and  III  show  that  these  subsidiary  conditions 
can  be  treated  in  minor  iterations  for  each  major  iter- 
ation step  of  the  natural  iteration  (NI)  technique. 

These  sections  also  present  computational  details  of 
calculating  the  second-order  transition  point  without 
losing  the  spirit  of  the  NI  method.  This  technique 


makes  most  use  of  the  superposition  expressions  which 
form  the  starting  point  of  the  NI  treatment. 

The  scalar  product  (SP)  expression  of  boundary  free 
energy”  has  lacked  a rigorous  proof.  Section  IV  pre- 
sents an  improvement  of  the  proof.  The  SP  expression 
is  combined  in  Sec.  V with  the  results  of  Secs.  11  and 
III  for  the  W-  and  double-square-approximations  to  cal- 
culate boundary  free  energies  o(10)  and  oul)  for  the  two- 
dimensional  Ising  model.  The  resulting  curves  in  Fig. 

3 leave  no  doubt  that  the  a values  calculated  by  the  SP 
expression  converge  to  the  exact  results  of  Onsager  and 
Fisher  Ferdinand  as  the  approximation  is  improved. 
This  convergence  supports  the  correctness  of  the  im- 
plicit assumptions  still  left  in  (he  analytic  proof  in  Sec. 
V,  and  now  guarantees  that  the  SP  expression  can  be 
safely  used  for  further  studies  of  the  boundary  struc- 
tures. 
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Structure  of  phase  boundaries* 
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The  boundary  between  plus  and  minus  spin  phases  of  the  (sing  model  is  studied  for  the  bcc  <110> 
boundary  The  iterative  calculation  due  to  Weeks  and  Gilmer  is  interpreted  as  a modification  of  the 
Natural  Iteration  computation  (with  subsidiary  conditions)  for  deriving  a free  energy  minimum.  A 
tetrahedron  is  used  as  the  basic  cluster  in  the  cluster-variation  scheme  for  this  problem.  Two  formulations 
are  presented:  one  uses  the  scalar-product  formulation  and  the  other  the  boundary  sum  method.  Previous 
results  of  an  order-disorder  type  phase  transition  within  the  boundary  are  confirmed.  Results  of 
calculations  on  the  boundary  excess  free  energy  and  the  boundary  profile  are  compared  for  the  pair  and  the 
tetrahedron  treatments. 


I.  INTRODUCTION 

In  a recent  publication.  Weeks  and  Gilmer1  proposed 
a novel  iterative  technique  for  numerically  solving  the 
density  profile  across  a phase  boundary.  They  pre- 
sented their  idea  based  on  the  original  Bethe2  method 
ol  treating  cooperative  systems  using  consistency  rela- 
tions. In  Sec.  II  of  the  present  paper,  we  point  out  that 
the  Weeks-Gilmer  (WG)  technique  fits  well  with  the 
cluster- variation  method  of  formulating  the  boundary 
structure  by  Cahn  and  the  present  author,3’4  and  that  the 
WG  iteration  technique  supplements  the  Natural  Itera- 
tion (NT)  technique  recently  proposed  by  the  present 
author. 5,3  The  iterative  technique  of  Sec.  n can  be 
generalized  to  larger  clusters.  In  Sec.  IV  we  calculate 
the  boundary  profile  of  the  bcc  (110)  boundary  using  a 
tetrahedron  as  the  basic  cluster. 

Along  with  the  boundary  profiles , the  excess  free  en- 
ergy of  the  boundary  is  calculated  in  Secs.  II  and  IV. 

In  order  to  compare  these  results  with  those  obtained 
by  the  scalar-product  formulation, 8,7  Sec.  HI  reports 
the  latter  using  the  tetrahedron  as  the  basic  cluster. 
Particular  emphasis  is  placed  on  the  phase  transition 
within  the  phase  boundary.  Examples  are  limited  to  the 
(110)  boundary  of  the  bcc  Ising  model. 

II.  ITERATIVE  COMPUTATION  OF  THE  PROFILE 
The  pair  approximation  formulation 

In  this  section  we  use  the  pair  approximation  of  the 
cluster- variation  method  to  calculate  the  density  (spin) 
profile  across  the  bcc  (110)  boundary  and  show  where 
the  WG  idea  fits  into  the  formulation  and  how  their  idea 
facilitates  numerical  computations.  We  work  with  the 
nearest-neighbor  interaction  as  in  WG  and  use  the  Ising 
model  terminology,  rather  than  ordered  alloys. 

The  crystal  structure  is  illustrated  in  Fig.  1.  We 
consider  the  plane  of  the  lattice  points  i-j—l  to  be  paral- 
lel to  the  boundary,  and  phase  I to  be  on  the  left  and 
phase  II  on  the  right.  A lattice  plane  parallel  to  the 
boundary  is  called  a '‘parallel’’  plane  for  short  and  is 
numbered  by  v,  as  in  Fig.  1. 

Two  kinds  of  nearest-neighbor  bonds  are  to  be  dis- 
tinguished, to  be  called  the  B bond  and  the  C bond  for 
short.  A B bond  is  within  a parallel  plane  and  examples 
are  i-l  andj-1  in  Fig.  1.  AC  bond,  such  as  i-m  in 
Fig.  1,  connects  two  lattice  points  on  adjacent  parallel 
planes.  A B bond  on  the  vth  plane  will  be  called  a Bv 
bond,  and  a C bond  connecting  points  on  the  vth plane  and 


the  (v+  l)th  plane  will  be  called  a Cp  bond.  The  number- 
ing system  p for  bonds  is  indicated  in  Fig.  1.  Actually, 
we  use  jx  = v in  numerical  computations,  but  it  is  better 
to  distinguish  them  in  the  presentation. 

Plus  and  minus  spins  are  denoted  by  j = 1 and  2,  respec- 
tively. The  distributions  of  spin  configurations  over  a lat- 
tice plane  and  over  bonds  are  defined  as  in  Table  I.  We  use 
distributions  of  nearest-neighbor  pairs  only  in  this  section. 

There  are  three  kinds  of  constraints  for  these  vari- 
ables: first,  the  set  of  normalizations 

HyBv.ii=Sj’c*.jm=l  for  all  i/’s  and  p’s;  (2.1) 

(.1  J.m 

second,  the  symmetry  relations 

^bv.ii  =yBv,n  for  all  v s;  (2.2) 

and  third,  the  consistency  relations 


:v,i = S y b Vtt  1 1 

i 

xv,i = ^2yBVtn 
l 

(2.  3a) 

Vti  ~ ^ yCUtJm* 

1 e *R  “ > y^  UeJm 

(2.  3b) 

m u 


In  Eqs.  (2.  3b),  remember  that  the  Cp  bond  connects 
lattice  points  between  the  vth  plane  and  the  (i>  + l)th 
plane,  as  shown  in  Table  I and  Fig.  1. 

The  cluster-variation  method  starts  with  the  free  en- 
ergy expression.  The  energy  of  the  system  is  written  as 


m-1  \ 

V 

FIG.  1.  Geometry  of  the  bcc  lattice  and  the  nomenclature  used 
In  defining  variables.  Thin  lines  are  nearest-neighbor  bonds, 
and  thick  lines  are  second-neighbor  bonds.  Section  II  uses  only 
thin  line  bonds. 
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TAHLK  I.  Definition  of  the  dial ribuiai Ion  variables.  The  sub-  where  the  £ operator  originates  Iront  the  Stirling  ap- 
srrtpts  u in  va„i(i  are  symmetric  as  seen  In  (2.2).  For  proxlntation  of  a factorial  and  is  defined  as 

y t i he  first  subscript  ( refers  to  the  left  side  of  the  bond 

as  shown  In  Fig.  1.  £(*)  * ■*  In*  — * . 


Configuration 
(-spin  on  the  I'th  plane 
i-i  on  a Bv  bond 

i-m  on  a Cp  bond, 
i on  the  rth  plane,  and 
hi  on  the  (c  > l)th  plane 


£- 2m£c„  + l^cu.ii] 


where  n is  the  number  of  lattice  points  on  each  of  the  parallel 
planes.  The  energy  parameter  for  the  lstng  model  is 

€,,*  + « when  is j. 


c,j  = -€  when  i=j. 


Distribution 

variable  In  (2.6),  JV  is  the  number  of  lattice  points  in  the  system 

and  2w  is  the  coordination  number.  For  an  inhomogene- 

ous  system  in  a bcc  lattice  (u>  = 4),  which  is  what  we  are 
Vav.ri  Interested  in,  we  modify  (2.6)  and  write  the  entropy  as 

♦ r££uvtJ)t££(*v.,.ji  (2.8) 

(2.4)  “‘lL'  “ J 

of  the  parallel  zL  [Siity.v.,,)-  i]-2g  j]  • 

g model  is 

The  four  x terms  in  (2.  8)  are  equal,  except  for  the  end  terms, 
but  are  written  separately  to  make  the  treatment  symmetric. 


(2.  5)  The  - I terms  in  (2.  8)  take  case  of  (u>  - 1)  in  (2.  6). 


In  (2.4),  the  parallel  planes  go  from  v=  1 through  p„ 
and  the  bonds  go  front  p = 1 through  p,  = v,  - 1. 

The  entropy  of  the  pair  approximation  of  the  cluster- 
variation  method  is  written  for  a homogeneous  phase  of 
S lattice  points  as*’* 


In  minimizing  the  free  energy,  we  want  to  treat  all 
y’a  as  independent.  For  that  purpose  we  need  Lagrange 
multipliers  for  the  following  relations: 

*►,<  ~ Evav.ll  = > 

I m 

**,i  = Z^c<u-i>.»i  = 5j.Vav.ii  • (2.9) 


1 -i  •'«,!  *W.'SV,I|.  v ' 

(2w- l)Xv  £(*,)- £(yu)+(w- I)  , (2.6)  * 1 

i d J We  write  the  Lagrange  terms  as 

_J 

r = J]  {E.tol[E(x ^av.„)-2  Z^yc«.i«]  + £va«v,i|2]Cycte-i>.»<  _ 52  (>’av,u  + }*".»)]  | 

-Y.  iE(alv.l  + “tv, I ~ aRv.(  ~ aRv.|)>'8v,ll  + ^52  (aR(v*n,«~  atv,i  Vc  wlw  / • 

( <1  I * 

In  the  first  line  of  this  expression,  yBv,„*y»v,u  .VcJ.jm  = US|X^,.*)7/'exp(-  0c,  J . (2. 13) 

guarantees  that  is  invariant  under  the  interchange 

of  i and  l.  ' Note  that  the  il  exchange  symmetry  of  in  (2.2)  is 

demonstrably  satisfied  by  the  expressions  in  (2. 12)  and 
The  free  energy  of  the  entire  system  F=  E-  TS  is  ^ 13), 

written  using  £ in  (2. 4)  and  S in  (2.  8)  and  adding  the 

Lagrange  terms  T (actually  - 2T)  in  (2. 10)  as  When  4>  in  (2. 11)  Is  a minimum  and  (2. 12)  and  (2. 13) 

/ \ hold,  we  can  show  that  <t>  reduces  to 

*■?“? -r*-2r+£>4-p-')  * a 

“•  . E/n  = 22^*o*  (2.14) 

+052icM-52yc»l2)  ■ (2-id 

' *•"  ' This  F is  the  free  energy  of  the  entire  system  including 

The  last  two  summations  are  the  Lagrange  terms  to  take  the  boundary  region, 
care  of  the  normalizations  (2. 1). 

The  Lagrange  multiplier  a’a  are  determined  from 
The  equilibrium  state  of  the  system  containing  the  ^ ^ Let  ug  de{lne 

boundary  is  derived  by  minimizing  * in  (2. 11).  In  dif- 
ferentiating with  respect  to  the  y’s,  we  regard  v’s  in  (0)  , . 

each  of  the  four  terms  in  (2.  8)  as  dependent  on  the  y’s  Sbm  ■2jv.Vav.ii  exP  ~ atu.i  + a«v.i  • 

written  explicitly  in  (2.3).  Minimization  leads  to 

Vav.ii  ■ .Va°v.ii  exp(i0*av"  atv.i  - atv.i  + anv,i  + aav.i).  sc*v,t  * 52ycltm  <*p(-  2aa(wl)f»)  , 

Vcu.ii" “ ycv.j"! exP(i^>cv  + 2atv,i  “ 2aa(v.,),»i)  > ^ 


SB*, , ■ 52  Van , , exp(-  aLV, , + aK„, , ) , 

I 

ScRv.l  ■£vc°wi«exP(-2a  aiv»i  >, »)  > 

m 

f>Civ,l  * Evc°('u-l).»l  exp(2at  (vl>,»)  • 


Vav..i  = <*v.<*v.i>1/8exP(-  0ei«>  « 


Substituting  these  in  (2. 9),  we  obtain 
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exp(i0X 

8»  - + aKi’.l  “ CXP(20*CM  ^ )ScR»t(  * 

exp(|0A  c(m-i>  - )Scil>l) 

■ exptj^Xj^-  + ®nt»ii)5Bk,i  • (2.16) 

The  novel  idea  of  Weeks  and  Gilmer1  is  to  eliminate 
the  x’s  from  (2. 16)  by  forming  ratios  of  i=  1 and  i = 2 
cases: 

e*Pj~  <W+  - exptta^jS^,,., 

exp(-aIK.2+a(B.,i)SB„l  exp(2atl,iJ)sCRl,tJ  ’ 

exp(-2ail>tt)SCI|>[|  _ exp(-  aL>- , + i)Sflw  i -7> 

exp(-  2 aR„% j)$c l*. z exP(—  aL»,i+  aR*.j'SB„, 2 

From  these  we  can  solve 

aRf,i~  ajj»,2=  (3ylilV - /lRt„)/8  , 

®i», i — a i M,t~  » (2. 18) 

where  the  A’ s are  defined  using  the  S’s  in  (2. 15)  as 

“ ln^cRi.,  I^Bm,  j/ Scrv.zSbv,  i)  > 

^it*,l|1(Sci»iA»,|/Scii’.A»li)  • (2.19) 

Since  the  y’s  satisfy  the  normalizations,  we  can  choose, 
for  example, 

«ik,2=  a*„,8  = 0 (2.20) 

without  loss  of  generality. 

The  iteration  procedure  is  divided  into  ‘major”  and 
“minor”  iterations.'  A major  Iteration  step  starts  with 
the  input  set  {x {alKll},  and  {<*„„_ Using  these, 
we  calculate  and  y££ jm  in  (2. 13)  and  determine  the 

output  set  {5}  by  the  minor  iteration  steps  to  be  de- 
scribed below.  When  the  output  a’s  are  determined, 
the  x’s  are  calculated  by  substituting  (2. 12)  into  the 
normalizations  (2. 1);  then  the  information  of  y,0”s,  x’s, 
and  a’s  leads  to  the  output  sets  xnd  {yC|ltim}  from 

(2. 12).  The  output  set  {£„,,}  is  obtained  from  (2.  9)  us- 
ing y’s.  Thus,  one  major  iteration  cycle  is  completed, 
1,1,1  {*►.(}>  M<l  {S„„(<}  are  used  as  the  input  set 

for  the  next  major  iteration  cycle. 

The  minor  iteration  step  is  as  follows.  The  input  for 
the  minor  iteration  are  the  sets  {y(0)}  and  {a}.  Using 
these  values,  we  calculate  the  S’s  in  (2. 15)  and  further 
A’s  in  (2. 19)  and  then  obtain  the  output  values  aR,tl  and 
aL».i  of  one  minor  iteration  cycle  from  (2. 18)  and  (2.20). 
This  output  set  a’s  is  then  used  as  the  input  on  the  right- 
hand  side  of  (2. 15)  for  the  S’s  in  the  next  minor  itera- 
tion cycle. 

We  can  prove  the  convergence  of  the  major  iteration 
by  showing  that  the  free  energy  function  <J>  in  (2. 11)  al- 
ways decreases  at  each  major  iteration  step.  The  ex- 
pression <1>  is  a function  of  the  variables  x’s,  y’s,  and 
a’s.  When  these  variables  take  the  output  values,  we 
write  <$  with  a caret,  and  we  let  <f>  denote  the  function  of 
the  input  values.  Then  using  a technique  similar  to  the 
one  used  in  Ref.  5,  and  assuming  that  the  constraints 
(2.  9)  are  satisfied,  we  can  derive 

* - 2 2Z  [ yBv,ti 

M7l 

+ ^ [ ycm.j"  ^n(ycu,lm/ycn'/iii)  *$Cn,/m  ~ ycu,lm\ 


+ 7 21  [*...<  tote... (/*...( ) + - K.i ] • (2.21) 

Each  term  in  the  square  brackets  is  nonnegative  due  to 
Gibbs’  lemma, 10  and  hence  <f>  always  decreases  at  each 
major  iteration  step  as  long  as  the  minor  iterations 
have  converged  and  the  continuity  relations  (2. 9)  are 
satisfied  by  both  the  Input  and  the  output  y’s. 

We  have  not  attempted  the  analytic  proof  of  the  con- 
vergence of  minor  iterations.  In  numerical  computa- 
tions they  do  converge,  however,  and  the  reason  for  this 
convergence  is  interpreted  as  due  to  the  averaging  na- 
ture of  (2. 15)- (2. 18). 

We  want  to  point  out  that  the  iteration  processes  of 
this  section  are  almost  exactly  the  same  as  the  itera- 
tion method  described  by  Weeks  and  Gilmer. 1 The  only 
difference  is  that  they  do  not  treat  the  major  and  minor 
iterations  as  rigidly  separated  as  in  this  section.  In 
the  method  proposed  in  this  section  the  output  set  {a} 
makes  the  continuity  requirements  of  (2.9)  satisfied  for 
the  given  {y<0>},  while  in  WG’s  method  the  output  {a} 
of  one  minor  iteration  cycle  is  regarded  as  the  output 
of  the  major  iteration  cycle  also.  The  fact  that  their 
iterations  converge  indicates  that  we  do  not  need  to  be 
so  strict  about  minor  iterations,  although  the  conver- 
gence of  the  latter  is  used  in  our  proof  of  4>  > $ in  (2.21). 

We  can  derive  equations  in  the  Kikuchi-Cahn  paper 
from  (2. 12)  and  (2. 13)  by  eliminating  the  x’s  and  a’s: 

3Wi  yc».n3lc(i»-i>.«  -/yKiV/i! 

yBv.Zt  y C(K- 1), 22^ Cm. 21  V*V,2/ 

g4a< =ycu,nycu,a/ycu,izycu,zi » 

=y'BM,  ll>'Bi/,22Alut  12  • (2.22) 

These  are  exactly  Eqs.  (2.8)-(2. 10)of  Ref.  3 for  the 
case  u>  = 4 and  u>,  = ui,  = 2 in  their  notation.  Although  this 
set  of  equations  is  equivalent  to  the  mathematical  prob- 
lem we  treated  in  this  paper  using  the  NI  method,  when 
written  in  the  form  of  (2.  22)  the  NI  technique  is  of  no 
help.  The  Kikuchi-Cahn  paper  solved  it  successively 
starting  from  one  end. 

It  may  also  ’.•«#  worth  pointing  out  that  the  present 
method  which  is  based  on  the  cluster-variation  formula- 
tion allows  us  to  calculate  the  free  energy  of  the  system 
in  (2. 14)  in  terms  of  the  Lagrange  multiplier  x’s  for 
the  normalizations.  This  is  in  contrast  to  the  WG  for- 
mulation which  is  based  on  the  original  Bethe  treatment 
and  does  not  lead  to  the  free  energy  value  in  easy  steps. 

The  formulation  in  this  section  was  used  in  calculating 
the  boundary  structure  for  the  interstitial-center  case 
and  the  atomic-center  case.  The  results  agree  with 
those  of  Kikuchi-Cahn,5  and  are  compared  with  those 
of  other  methods  in  Sec.  IV.  3. 

III.  TETRAHEDRON  TREATMENT  OF  BOUNDARY 
FREE  ENERGY  USING  THE  SCALAR  PRODUCT 
FORMULA 

The  iterative  calculation  of  the  boundary  profile  pro- 
posed by  Weeks  and  Gilmer1  and  reinterpreted  in  the 
previous  section  in  terms  of  the  NI  calculation  suggests 
that  a similar  technique  can  be  used  with  a larger  clus- 
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ter  than  the  patr.  This  is  done  In  Sec.  IV  using  a tetra- 
hedron as  the  basic  cluster  for  the  bcc  (110)  boundary. 


As  preparation  for  the  formulation  in  Sec.  IV,  we 
solve  the  bulk  properties  of  the  bcc  istng  lattice  using 
the  tetrahedron  approximation  in  Sec.  ill.  The  bulk 
calculation  is  reported  here  as  an  introduction  to  the 
more  Involved  formulation  for  the  boundary  in  Sec,  IV, 
and  also  because  it  leads  to  the  scalar-product  (SB) 
evaluation  of  the  excess  boundary  free  energy  o.  The 
S1J  formulation  is  in  Sec.  ill.  2 and  the  results  are  com- 
pared with  those  of  Sec.  IV  in  Sec.  IV.  3. 


The  tetrahedron  approximation  for  the  bulk  reported 
here  has  been  used  before,  some  results  were  pub- 
lished, and  the  method  was  briefly  sketched  previ- 
ously, ” There  are  also  three  papers  published  by  a 
Russian  group  using  the  same  cluster. 14 


1.  Bulk  phase  using  tetrahedron 

The  tetrahedron  we  use  in  this  calculation  is  i-j-k-l 
in  Fig.  1,  which  is  not  a regular  tetrahedron  but  is  made 
of  four  nearest-neighbor  bonds  i-k.  i-l,  j-k,  i-l  and 
two  second- neighbor  bonds  i-j  and  k-l.  The  distribu- 
tion variable  for  a lattice  point  is  written  as  and  that 
for  a nearest -neighbor  pair  is  written  as  as  can  be 
understood  from  Table  I.  Besides  v and  v,  we  need 
three  other  kinds  of  distribution  variables:  for  a 

second- neighbor  pair  (like  i-j  in  Fig.  1),  utn  for  a tri- 
angle of  the  shape  i-j- 1 in  Fig.  t,  and  rr(M1  for  a tetra- 
hedron i-j-k-l  in  Fig.  1. 


The  entropy  of  the  system  is  written  as 


(3.  1) 


S-*lnO  , 

where  s)  is  the  number  of  ways  in  which  a system  can  be 
constructed  for  spec  if  i u ,et  of  distribution  variables 
{N'n*M  ftnd  is  written  for  a system  of  .V  lattice  points  as 

0 {Triangle  {point  {Tetrahedron 
x {Second-neighbor  pair  z0}'5 
x {Nearest-neighbor  pair  y,,};4  . (3.2) 


The  meaning  of  the  curly  bracket  notation  is  0X1)181111x1 
in  any  of  the  cluster-variation  publications  previously 
referred  to.  *•* 


When  the  expression  for  n is  given  as  in  (3.2),  it  is  a 
simple  procedure  to  write  the  free  energy  of  the  system 
in  terms  of  the  distribution  variables  and  then  minimize 
it  to  obtain  the  equilibrium  state.  In  the  NI  procedure 
we  treat  the  h(>„’s  as  independent  and  all  other  vari- 
ables as  dependent.  Following  steps  similar  to  those 
presented  in  Sec.  Ill  of  Ref.  5 to  keep  the  symmetry  of 
the  formulation  explicit,  we  arrive  at  the  "superposi- 
tion expression 


'U»l  • " 0*1 


(3.3a) 


M'i  m B *’ 

where 

" nil  ■ e*P<-  d«(/,|K«(HMl>,M,|,»,|/)1/,(vl  v<v,AI)1/H 

* ^ V| i.v^ «.v«|.vay 1 /44  , (3.3b) 

The  energy  parameter  is  defined  as  the  energy  per 


(3.4) 


tetrahedron  and  is  written  m terms  of  the  nearest- 
neighbor  energy  parameters  tu  of  (2.  5)  as 

*u*i  * (*u  * *ii  * *«i  4 **/'  ® • 

The  definition  of  in  this  |»apor  is  different  from  that 
in  Eq.  (9)  of  Ref.  13;  in  the  latter,  is  defined  as 
the  energy  contained  in  a tetrahedron  and  the  factor  of 
6 results  because  a nearest -neighbor  bond  is  shared  by 
six  tetrahedra.  The  present  definition  is  preferred  in 
anticipation  of  the  formulation1*  in  which  the  second- 
neighbor  interaction  energy  is  to  be  taken  into  account. 


The  Lagrange  multiplier  x in  (3.  3a)  is  for  the  normal- 
ization condition 

T u 


i 


(3.  5) 


and  as  is  usually  the  case,5  is  equal  to  the  free  energy 
of  the  system,  per  lattice  point: 


F A'*X  , 

when  the  free  energy  is  a minimum. 


(3.6) 


Because  of  the  symmetry  of  the  bulk  phase  which  is 
contained  in  the  rcu,,  expression  in  (3.  3),  we  do  not  need 
additional  symmetry  constraints  nor  the  Lagrange  mul- 
tiplier o's  as  are  used  in  Sec.  11  and  in  Ref.  6.  The 
Natural  Iteration  procedure  works  when  we  combine 
(3.  3)  with  the  geometrical  reduction  relations: 


UU»  = £nO*m 

i « 


(3.7) 


For  the  measure  of  convergence  of  the  Nl  procedure, 
we  take 


£ |lmc<%1,-lnH>>,| 


(3.8) 


where  the  superscript  («)  denotes  the  nth  iteration  step. 
This  test  value  A*"'  decreases  below  10'5  when  n is 
about  400  iterations  for  almost  any  temperature  in  the 
ordered  phase  in  this  particular  problem.  The  differ- 
ence of  logarithms  in  the  measure  of  the  test  as  in  (3.8) 
is  preferred  to  the  straight  difference  of  re's,  because 
the  former  yields  same  number  of  digits  accuracy  for 
all  re’s,  which  is  the  accuracy  we  need  in  the  scalar- 
product  evaluation  of  the  boundary  free  energy  as  was 
discussed  in  Ref.  6. 


2.  Scalar  product  evaluation  of  o 

The  excess  free  energy  a of  the  (110’'  boundary  can 
be  calculated  using  the  bulk  phase  results  of  Sec.  111.  1 
and  the  SB  expression.  A (110)  plane  of  the  bcc  lattice 
is  drawn  in  Fig.  2.  The  lattice  constant  a is  indicated. 


The  SB  expression  for  the  boundary  excess  free  en- 
ergy has  previously  been  derived  and  proved. ls’1*  Ap- 
plied to  the  present  problem,  we  write 


exp(-m>(M0)do*  Y2)«  £ OR..! I f>, {«,,,} /'n{«,„}|,/‘, 

l“4'*'  (3.9) 


rhere  n*  \2  is  the  area  per  atom  in  the  "parallel 
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bonds,  and  AB  and  i'l>  as  the  second-neighbor  bonds, 
as  Indicated  in  Fig.  2 of  this  paper.  The  result  is 

1){m,„}  {Second -neighbor  pair 

x {Nearest-neighbor  pair 

^{Triangle  n((,}^{Polnt  ij;1  , (3.10) 

where  the  curly  bracket  notation  has  the  same  meaning 
as  that  m (3.  2)  or  as  tit  those  in  Refs.  8 and  9. 

Ckir  experiences  in  Rets.  7,  17.  and  6 tell  us  that  we 
can  use  the  schematic  Information  in  t)  of  (3.  10)  when 
we  write  the  probability  function  appearing  in 

(3.9).  It  is  written  as 


/'ll"..*}  i/*f..fi„ pi!, />u,y 


(3. 11a) 


KIU.  -.  4\llo  pl*iw  *>l  ihf  Kv  lattice.  I* he  iwaivai-iun^hbor 
•H»nd  ami  th«'  Mi  dml  m inhlKir  boiwt  aiv  drawn  with  a thin  line 
.tml  a ihii’k  hiw.  itM|Hvi iveh 


plane,  n is  the  total  number  of  lattice  points  in  the 
plane,  and  o(u#(  is  the  excess  boundary  free  energy  per 
unit  area.  On  the  right-hand  side,  />i{m,„}  is  the  prob- 
ability t ha',  a nitt  plane  specified  by  the  set  of  vari- 
ables {»,,,}  is  realized  in  phase  1,  and  il{n(M}  is  the 
number  of  ways  the  plane  can  be  constructed  when  the 
set  {«, „!  is  specified. 

The  weight  (actor  l){ii,y,}  can  be  calculated  by  slightly 
modifying  the  scheme  ti sed  in  See.  G of  Ref.  8.  We 
can  follow  the  steps  of  Eq.  (G.  1)  if  we  regard  bonds  AC 
and  .-if)  in  Fig.  8 of  Ref.  8 as  the  nearest -neighbor 


/•i.."  II  <«/)*’>  ••«!/*  . 

t.i 

pi.**  n / «»* 

4,»  V '/*  >**  ' it  • 

/>!.«•  n (*ia>)»* . 


(3.  lib) 


The  double  asterisk  is  the  Fortran  notation  for  "raised 
to  the  power  of,  and  is  used  here  to  avoid  subscripts  of 
a superscript  in  typesetting.  The  quantities  with  the 
superscript  (1)  in  (3. 11),  such  as  «/},',  are  those  for  the 
bulk  (1)  phase.  Note  that  the  powers  on  the  />’s  on  the 
right-hand  side  of  (3. 11a)  correspond  to  those  in  (3.  10). 

When  we  use  (3.  10)  and  (3. 11)  in  (3.9),  we  derive 


o<ll0)/ta*/vT=2  E £(*,,,)- E£(V-E*(.v,»)-E  *(>„)♦  Y>Iv,kEjC(v,)1 

i»j»*  i »i  it*  .f,fe  H i i t J 

- 2 E "ut  * iE  lhijy + E in'(*  • E >/» 

<•>.»  i.i  i,i  rX 

- gfE  'i  *E  ► E '* * dx  ( 1 - E > 2 E (at,  + Oy,  - o«(  - »*>)«..,  . (3.12) 

L 1 > * J ' 1,4.1  / 1,1,* 


The  £ operator  is  defined  in  (2.  7).  The  quantity  with  a 
caret  is  defined  as  the  geometric  average  of  the  corre- 
sponding bulk-phase  quantities  for  phases  1 and  II  as 

*iz  ■ U,*{‘*/**’>,/*  . etc.  (3.13) 

The  Lagrange  multipliers  a„  are  Introduced  to  take 
care  of  the  symmetry  of  Vy,: 


.v*-£".!.*E"<*,.  (3.14) 

In  minimizing  o in  (3. 12),  we  use  the  geometric  rela- 
tions 

*,1-E-.1*.  *|-E“i/»  (3.15) 

• 4.* 

together  with  (3. 14).  The  minimization  leads  to 


*hi*  = **ii*  exp(l$  x + a„  + Oy,-  a„  - «*y)  . (3. 16a) 


to.  n ^liVt*.',*)1'* 


(3.  16b) 


When  (3. 16)  holds  and  o in  (3. 12)  is  a minimum,  we  can 
show  that 

x *0<iio>",,v'2  . (3.17) 

which  is  the  excess  boundary  free  energy  per  lattice 
point  in  the  plane,  it  is  easy  to  see  that  v - 0 when  the 
bulk  phase  is  in  the  disordered  state  so  that  m",’  m,1}}’ 


The  Lagrange  multipliers,  the  o(,’s,  are  determined 
from  the  constraints  in  (3. 14).  After  transformations 
similar  to  those  in  (2. 15M2. 19),  we  arrive  at 
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,,,  i “m  * “mexP^~ 

.Jv.iSsKi;/  • 


• <*)•  - 


®n  ^ *** 3 °*  *iis“4ji*  (3.19) 

We  can  solve  the  set  o{  simultaneous  equations  by  the 
Natural  Iteration  technique.  We  need  the  major  and 
minor  iterations  as  was  the  case  in  Sec.  11.  However, 
since  the  constraint  equation  in  (3. 18)  in  a fourth  or- 
der algebraic  equation  for  expiA,,),  the  minor  itera- 
tions can  be  done  by  the  Newton- Raphson  iteration 
method  without  trouble. 

3.  Phase  transition  within  the  boundary 

As  was  reported  in  Ref.  17,  the  scalar-product  treat- 
ment of  the  excess  boundary  free  energy  c(no)  using  the 
pair  approximation  leads  to  a second-order  phase  change 
in  °<uo>-  We  niav  c*l* this  transition  temperature  T„. 

We  expect  a similar  phase  change  in  the  tetrahedron 
(bulkMrlangle  (boundary)  treatment  of  Sec.  111.2. 

In  the  triangle  formulation  of  Sec.  HI.  2,  we  see  that 
there  are  two  long-range  order  parameters: 

( , ■ x,  - jr|t  (aaMU|-Mul.  (3.20) 

The  transition  point  T,  can  be  determined  as  the  point 
at  which  the  following  2x2  determinant  vanishes  in  the 
disordered  phase: 

det||8,o(uol/8t,8«J||  = 0 . (3.21) 

The  determinant  can  be  worked  out  explicitly  and  the 
equation  takes  the  form 

k-±  * (l-i-±io. 

V y 1 1 fi|i/  \miii  wmi/ l \ .'n  2*ii/J 

(3.  22) 

It  may  be  commented  that  in  the  disordered  phase  A„  in 
(3. 19)  vanishes  and  we  do  not  need  the  minor  iteration 
of  (3. 18).  When  we  solve  the  disordered  phase  In  Sec. 

III.  2 and  evaluate  Eq.  (3.22),  we  do  arrive  at  the  sec- 
ond-order transition  point  T ,,  qualitatively  confirming 
the  previous  conclusions  of  the  pair  approximation. 1T 
The  results  are  discussed  in  Sec.  IV.  3. 

IV.  TETRAHEDRON  TREATMENT  OF  BOUNDARY 
FREE  ENERGY  USING  THE  SUM  METHOD 

1.  Free  energy  and  its  minimization 

We  are  now  in  a position  to  continue  from  Sec.  II  and 
extend  the  pair  treatment  of  the  boundary  sum  into  the 
tetrahedron  treatment.  From  Fig.  1,  we  see  that  there 
are  three  kinds  of  tetrahedra  in  an  Inhomogeneous  sys- 
tem: one  pointing  towards  the  left  (i-j-k-l),  one  point- 
ing towards  the  right  U-j-l-ni ),  and  a third  one  (i— »»— 
/-«i).  We  write  the  distribution  variables  correspond- 
ing to  these  three  as  »**„*,,  fvKrtWm,  and  re- 

spectively, where  v and  p indicate  the  position  as  seen 
in  Fig.  1.  Of  the  four  subscripts  such  as  ijkl,  the  first 
two  and  the  last  two  connect  the  second-neighbor  lattice 


points.  In  each  of  the  second-neighbor  subscript  pairs, 
the  first  one  is  either  above  or  on  the  left  side  of  the 
pair.  The  ie’s  satisfy  the  following  symmetry  relations 
among  subscripts: 

“'h-,i«i  “ M’i*.)i»i  • 

M’S», Ilia  > 

M>cu,iiii« “ 8'cm, i«u  • (4.1) 

As  subclusters  of  these  tetrahedra,  we  see  there  are 
five  kinds  of  triangles,  for  which  the  distribution  vari- 
ables are  written  as  «s«.i,n  »nfl 

“l  »,!»«■  These  five  triangles  can  be  identified  by  follow- 
ing the  subscripts  and  the  corresponding  points  in  Fig. 

1.  Among  the  triplet  subscripts  like  ijl,  the  first  two 
connect  a pair  of  second- neighbor  points.  The  n’s  are 
reduced  from  w’a  by  the  following  two  sets  of  geometric 
relat  ions: 


» "i 

MSu,li»l  = 22k’c».Ii<I»  * ^M'sv,UI«  • 
n i 

ul>.lin*iCl,'ct.l<lK'^,*'UnlUila  • (4.3) 

I / 

The  symmetry  relations  among  the  re’s  in  (4. 1)  are  in- 
herited by  m's  through  the  geometric  relations  (4.2)  and 
(4.  3).  The  equations  among  re's  in  (4.  3)  necessitate 
Lagrange  multipliers,  as  shown  below. 

Two  kinds  of  distribution  variables  for  a second- 
neighbor  pair  are  written  as  (which  is  within  the 

eth  plane)  and  which  connects  / of  the  nth  plane 

and  in  on  the  (v+  l)th  plane.  Together  with  n\  n,  and  t 
we  also  use  v and  y defined  in  Table  I. 

Since  we  use  si’s  as  the  Independent  variables,  the 
energy  of  the  system  in  (2. 4)  is  rewritten  in  terms  of 
re's  as 

r **  “«  -i 

£ = 2n  ^ JIM’n-,ii»i + X-  M’a«,u»i  + <C*’cu.u»i|  • 

l,J, ft,  I LiS  *M  »•!  J 


where  cOII,  the  energy  per  tetrahedron,  is  the  same  as 
the  expression  (3.  4)  for  the  bulk  phase.  The  multipli- 
cative factor  ii  in  (4.  4)  is  the  number  of  lattice  points 
within  a plane  parallel  to  the  boundary;  it  should  not  be 
confused  with  the  name  of  a lattice  point  n on  the  (e  + l)th 
plane  in  Fig.  1. 

The  0 expression  for  this  case  is  derived  by  modify- 
ing 0 in  (3. 2).  For  example,  the  12  triangles  appearing 
in  (3.2)  are  not  all  the  same  in  the  present  case.  By 
counting  the  frequencies  of  occurrence,  we  make  the 
following  replacement: 

{»‘lllVv  ~ {“Bi>,II|}*{mI*>>,I/*}i|{mPi>.IJ«i}« 

*{m*»,iii»}i!{mi »,»»■}»  • (4.5) 
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Applying  similar  replacements  to  other  factors,  we  ar- 
rive at 


Sblv.iii  * Y1w(li!,uh  exp(-  aLu-i).m  ~ at(»-i),*iil  > 


0=11  [{“B'.utii**  *»,<**}« 

vorit 

* {ul»,(nmt*  {**,<}«)  [ tjtiiZ  {wRy, </;«)*  {tt'c*.l’il»)» 

(4.6) 

The  free  energy  of  the  entire  system  can  be  derived 
using  (4. 4),  (4. 6),  and  S=  k lnfl  as  in  (3. 1).  We  will 
skip  writing  the  free  energy  expression  to  save  space, 
and  go  directly  to  the  relations  obtained  by  minimizing 
the  free  energy  with  respect  to  w^>u,m,  and 

“'on, i««.: 

Wtr,Hkl  =U>L*.I1H  exPH0*ii.  + aBiMJI  ~ al(it-l),tlJ 
- aK»-l>,»ll]  > 

WKr,Ulm  = wRv,lJlm  exP[l0*Ry  + aSu,lml  + aRu.lmJ  ~ aBR.IJli  » 
WC».lnlm  " “'cu.l’ll*  eXP(2  0*Cl»  + aLu,  I urn  * °L  >,ln  “ aKu.lml 

” ®ai*ii*il  > (4.7) 

where 

“’tv.i/ai  = expi-  ^*U»|)(“a».Ol“»i>.IM“t(»-l).»HuL(i.-l)^u)1/J 
X (*^|Xi>,yX»,|X^.lt,)  (^ai>,lizC(M*l>,»l^"I/< 

x ^>'a*,n>’al-J^yc(»-^^.»lJ'c(l«-l).»r)'I/,  > 

WRv,IJ)m  - exp(-  Ptf)lm)(uBi’,IJIur‘v,IJmullu.lmluJlu.lmj)1/* 

X ^f.l^M^v.l^K.1,*)  *^zBv,UzC 

* (ya». nyai’.iiycu.i.j'cu,/-) 1/6  > 

“’c  u.  I"l*  = exP(~  3«|nl*)(«B  u,  lnluR  li.  1*1  “t  m,l*muL  «,  l*n  )'/2 

X (xVliX„t  |x1»i,ii*'i*i,„)1/,4(zCm>  fmZc  mimY1** 
^^yBf.liy  Bt¥*lt,miC>Clt,lmyCu.liXllt  • (4.8) 

It  can  be  seen  that  the  three  mj,0”s  in  (4.  8)  reduce  to 
«i<0’  for  the  bulk  phase  in  (3.  3b)  when  the  system  is 
homogeneous. 

The  Lagrange  multipliers  \L„  \R„  and  XCu  in  (4.  7) 
take  care  of  the  normalizations: 

wL»,lll)l-  wRr,iilm " u,C».t*lM=  1 • 

(4.9) 

Similar  to  (2. 11)  and  (2. 14)  for  pair  treatment,  in  the 
present  case  the  free  energy  of  the  entire  system  in- 
cluding the  boundary  is  written  as  the  sum 

- *"  li”  +tJ  ^Hv+  Cm  • (4.10) 

* I'-S  l»M  U*1 

When  we  compare  this  equation  to  (3. 6)  for  the  bulk,  we 
see  that  each  \L„  \R„,  and  \Cu  reduces  to  x/3  when  the 
boundary  is  removed.  That  is  the  reason  why  the  0\L, 
term  in  (4.  7)  is  divided  by  2 in  contrast  to  the  factor 
J-  in  (3.  3a). 

The  Lagrange  multipliers  a’s  in  (4.  7)  are  determined 
by  the  three  continuity  relations  in  (4.  3).  It  is  conve- 
nient to  define 


Sbrv.iii  * ^“'».oi*exp(aR  U, Ini  + » 

^RRy,imi  * '^*WRvtiJlm  exP[  &R u,  imj  ~~  t 

SrC  u,  l mi  * ^^’cuelnliH  WPlfllLli,|«»i  + aL  u,lmn  ” aRut  ,nll  » 
■ 

^tCd.ln*5  S WCu,lnlm  eXP(ffli»,l*n  “ aRu,l*l  ~ a R u.  I nil  I 

* H—lWL  0<*I).y*(il  eXP[  aB  <y+l  ),imn  ~ aL  »,(*/]  * 

(4.11) 

Then  we  can  write  the  first  equation  in  (4.  3)  as 

~ exp(2^li>  + aBR,HI^BLi’.IJI 

= exp(jSXRv  — otB„t1jt)SBRi,t{ji  . (4.12) 

In  treating  this  set  of  equations,  we  can  follow  Weeks 
and  Gilmer’s  idea*  and  eliminate  x’s  by  forming  ratios, 
as  was  done  in  Sec.  II.  Further  we  note  that,  since  the 
w’b  satisfy  the  normalization  relations,  one  of  the  eight 
equations  (since  each  of  i,  j,  and  l takes  two  values  1 
and  2)  in  (4. 12)  is  redundant.  This  redundancy  also  oc- 
curs in  the  rest  of  the  equations  in  (4.  3).  Thus,  we 
can  choose  without  loss  of  generality,  for  example, 

aBv.  222  - *Rn,222  = a£ii,222  = 0 . (4.13) 

By  forming  ratios  from  (4. 12)  and  using  (4. 13),  we 
obtain  the  following  set  of  equations  to  determine  the 
a's: 

aB».U I = 2 ln(‘SflR.., (/I'Sg z&fSgRv, 222^8 />,(/[)  , 

aR  mi*i  = 5 ln(SRC u> ,mj SRRVi  zzz/Srcu  ,22zSrRv,  Ini  ) > 

ata,lfi*  = 2 222/S£I,  (b.d,222Si,Cb,(b„)  . 

(4.  14) 

From  (4.11),  we  see  that  SBL*,ul  and  SBR^UI  are  in- 
variant under  the  exchange  of  i and  j,  and  thus  as„t(iI  in 
(4. 14)  retains  this  symmetry.  Because  of  this  sym- 
metry in  aBy,in,  when  the  input  mi’s  satisfy  the  sym- 
metry properties  in  (4. 1),  the  symmetry  is  guaranteed 
to  be  inherited  by  the  output  mi’s. 

2.  The  natural  iteration  computations 

We  solve  the  simultaneous  equations  in  Sec.  IV.  1 us- 
ing the  Natural  Iteration  method.  As  was  mentioned 
following  (2. 20)  in  Sec.  II,  the  iterations  are  divided 
into  the  major  and  the  minor  ones.  A major  iteration 
starts  with  the  input  data  set  and 

{<**>.»  aRu>aL  »}•  We  calculate  w(£\  ui£0>u}  from 
(4. 8).  Then  we  do  the  minor  iterations  to  satisfy  the 
subsidiary  conditions  (4. 14).  After  the  minor  itera- 
tions converge,  the  output  is  written  as  {6B„,  fi„B,  Siu}. 
We  use  these  together  with  the  m <0”s  in  (4.8),  which 
have  already  been  calculated,  in  (4.  7)  to  derive  the  out- 
put set  {®tB,  m)^,,  *Cu}  together  with  the  set  {xtB,  xRb, 

\c  The  output  w’b  lead  to  the  next  input  set  {ms„,  . . . , 
“li*}  thus  closing  one  major  iteration  cycle. 

In  Cahn-Kikuchl’s  work*14  it  was  discovered  that  two 


1 
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cases  should  be  distinguished;  one  is  the  interstitial- 
center  (IC)  boundary  in  which  the  center  of  the  boundary 
is  between  two  parallel  planes,  and  the  other  is  the 
atomic-center  (AC)  boundary  in  which  the  center  of  the 
boundary  is  on  one  of  the  (>arallel  planes.  We  did  cal- 
culations for  these  two  cases.  In  the  1C  case  we  took 
the  number  of  parallel  planes  r,  20  (the  center  being 
between  i - 10  and  11),  and  in  the  AC  case  we  took  e, 

= 19  (the  center  at  r 10). 

In  both  cases  the  initial  condition  for  the  iteration 
was  chosen  such  that  the  left-hand  side  of  the  boundary 
center  is  in  the  bulk  i phase  and  the  right-haiui  side  in 
the  bulk  ll  phase.  When  some  variable  is  right  at  the 
center,  we  used  the  averaged  value  between  the  two  bulk 
phase  values.  The  initial  values  of  o’s  are  chosen  as 
zero  all  through  the  system. 


(N 


C 


In  carrying  out  this  iteration  calculation,  several 
comments  are  in  order.  These  comments  are  important 
when  a system  is  as  large  as  the  present  one,  in  which 
the  total  number  of  independent  i< ’s  and  or’s  is  about  700. 
(Actually  a system  much  larger  than  the  present  one  has 
been  successfully  calculated  bv  Sanchez  at  UCI.A1*  using 
the  Natural  Iteration  technique. ) The  first  comment 
concerns  the  minor  iteration  procedure.  Equation  (4.14) 
corresponds  to  the  minor  iteration  equation  (2. 18)  in  the 
pair  case  and  to  Eq.  (2.20a)  in  Ref.  6.  In  Sec.  II,  Eqs. 
(2. 17)- (2.  19)  work  and  converge  nicely  w ithout  trouble. 

In  the  present  case,  however,  when  we  start  from  the 
initial  set  a-0,  the  output  5 of  the  iterations  from 
(4.  14)  diverges.  This  was  not  completely  unexpected 
because  the  convergence  of  the  minor  iteration  is  not 
guaranteed,  although  the  convergence  of  the  major  itera- 
tion can  be  analytically  proved  as  in  Sec.  II.  In  the 
present  case  the  convergence  was  achieved  by  the  follow- 
ing trick.  Instead  of  using  the  output  & from  (4.  14) 


MAJOR  ITERATIONS,  n 

FIG.  3.  The  test  mraaure  A'"’  ill  (4.  It!)  of  the  NI  calculation  of 
ihe  bee  \1 10 ' boundary  sum  method  plotted  against  the  major 
iteration  n in  a semilog  scale. 


MAJOR  ITERATIONS  n 

FKi  4.  The  excess  five  energy  and  the  extrapolated  value 
using  (4.  18)  calculated  using  the  NI  method. 


directly  as  the  input  to  the  next  minor  iteration  cycle, 
we  use  the  average  as  the  next  input: 

l»PUt  “ + ^output  ) . l*l.  15) 

This  brings  the  next  input  closer  to  the  previous  input 
than  using  <Soutput.  and  it  works.  Depending  on  the  cir- 
cumstances, we  can  conceive  of  many  modifications  of 
the  scheme  in  (4. 15).  In  the  work  reported  in  this  pa- 
per, the  number  of  minor  iterations  at  each  major  iter- 
ation step  was  about  60. 


The  next  important  consideration  we  can  make  use  of 
is  the  extrapolation  of  the  major  iteration  steps.  For 
the  measure  of  convergence  of  the  major  iteration  we 
used  the  following  quantity: 


A'"’=  L £ {Ih'KYm, 

iJik.  I v or  u 


4.  I ...inl  I 

+ I '*  Rp.l|*l  N«»,  (J*ll 


(4.  16) 

Different  from  (3.8),  we  did  not  take  the  logarithms  of 
re’s  in  (4. 16),  but  it  is  adequate  for  our  present  pur- 
poses. This  test  value  a'"'  decreases  logarithmically 
as  illustrated  in  Fig.  3 for  kT/e-  3.0,  and  we  can  prove 
that  A,n'  is  linear  in  the  semilog  scale  for  large  num- 
ber ii  of  the  major  iteration. 


In  the  same  sense,  a also  approaches  the  converged 
value  o'”’  exponentially  as  is  seen  in  an  example  of 
kT/t~  3.0  in  Fig.  4.  Using  this  property,  we  can  ex- 
trapolate o starting  with  the  assumption 


(4. 17) 


where  c*"’  is  the  value  of  the  boundary  excess  free  en- 
ergy at  the  nth  major  iteration.  It  is  a simple  matter 
to  derive 

(a^T/U^'"-'"’- a„o'"')  , 

where 

> (4.18) 

The  extrapolated  value  o'”’  is  a function  of  n and  in  since 
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the  exponential  decay  relation  (4. 17)  does  not  hold  ex- 
actly unless  n is  very  large.  As  we  see  in  Fig.  3,  A*"’ 
is  still  slightly  curved  in  the  semilog  plot  at  n-  300. 

The  extrapolated  value  a*”’  is  also  plotted  in  Fig.  4 as 
a function  of  n.  ft  is  almost  flat.  This  plot  is  based 
on  the  choice  of  m = 5,  but  other  values  (in  = 10  and  20) 
give  practically  the  same  extrapolated  value  for  o'“\ 

The  same  extrapolation  idea  was  used  in  the  estimate 
of  the  profile.  Figure  5 shows  how  x¥tl  for  t/=  11  and  12 
are  extrapolated  using  the  same  scheme  as  (4. 18). 

Actually,  a much  more  desirable  scheme  is  to  ex- 
trapolate the  entire  set  t0 

n — <*>  starting  from  a finite  iteration  step  n using  a for- 
mula similar  to  (4.  18).  Such  an  extrapolation  has  been 
tried  and  found  successful  in  other  applications  of  the 
Natural  Iteration  method,  but  for  no  apparent  reason 
has  not  vet  been  applied  to  the  present  problem. 

3.  Comparison  of  results 

Figure  6 compares  the  results  of  the  excess  free  en- 
ergy a of  this  paper  and  previous  ones.3, 17  Of  the  four 
curves  of  o,  the  upper  one  indicated  as  PAIR  S. -P.  is 
the  result  of  the  scalar-product  formula  using  the  pair 
approximation  in  Ref.  17.  It  has  the  second-order  phase 
transition  at  T,  of  Table  II,  as  shown  by  a cross  in 
Fig.  6. 

The  broken  curve  named  PAIR  SUM  is  the  result  of 
Sec.  11  of  this  paper  and  is  the  same  as  the  result  shown 
in  the  curve  (2),  Fig.  9 of  Ref.  3.  This  curve  is  for  the 
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MAJOR  ITERATIONS,  n 

FIG.  5.  The  plus  spin  density  at  v-  11  and  12  planes  and  their 
extrapolated  values. 
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kT/e 

FIG,  6.  The  boundary  excess  free  energy  a for  the  bcc  <110) 
boundary  calculated  by  four  different  methods.  See  the  text. 

IC  boundary  and  is  lower  than  the  AC  value.  The  differ- 
ence of  the  two  boundaries,  dAC-ic>  is  plotted  at  the  bot- 
tom of  Fig.  6 in  a broken  curve. 

The  solid  curve  marked  TRH  S.  -P.  is  the  result  of  the 
scalar-product  method  in  Secs.  III.  2 and  III.  3 of  the 
present  paper.  The  second-order  phase  transition  point 
is  again  indicated  by  a cross  at  the  value  of  T„  of  Ta- 
ble II. 

The  dot-dash  curve  which  is  the  lowest  and  marked 
TRH  SUM  is  the  result  of  the  method  in  Sec.  IV  of  this 
paper,  calculated  for  the  IC  boundary.  The  AC  case 
was  also  calculated  and  the  difference  aAC-oIC  is  again 
plotted  by  a dot-dash  curve  at  the  bottom. 

The  first  general  conclusion  we  can  derive  from  the 
comparison  of  these  four  curves  in  Fig.  6 is  that  they 
agree  well  qualitatively  and  thus  guarantee  that  our 
methods  are  sound.  The  second  conclusion  we  derive 
is  that  the  sum  method  of  Secs.  II  and  IV  gives  lower  o 
values  than  the  corresponding  approximation  of  the 
scalar-product  results,  in  Ref.  17  and  Sec.  IV.  The 
intrinsic  reason  for  this  is  not  known  at  this  writing. 


TABLE  II.  Transition  point  data.  The  value  with  an  asterisk 
was  previously  reported  in  Table  I of  Ref.  13. 


Method 

Curie  point 
in  the  bulk 
kT,/t 

Transition 
point  in  <r 
*T„/€ 

<T(a,/V2)/« 
at  T „ 

Pair,  Ref.  17 

6.9521 

2. 8854 

3.  5992 

Tetrahedron, 

this  paper 

6.4907* 

2.4918 

3.7267 

i 
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LATTICE  PLANES  v ACROSS  THE  BOUNDARY 

FIG.  7.  Profiles  of  spin  density  for  the  bcc  <110)  boundary 
calculated  by  the  pair  method  (dotted  curve)  and  by  the  tetra- 
hedron method  (solid  curve). 


z 

* 
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FIG.  H.  Thickness  /,  of  the  boundary  defined  by  (4.  19)  for  the 
bcc  (110)  boundary  calculated  by  the  two  methods. 


The  profiles  of  the  IC  boundaries  are  plotted  in  Fig. 

7.  They  arc  for  temperatures  t-T/c  = 2. 0,  3.0,  4.0, 
and  5.0  from  the  sharpest  to  the  most  gradual  in  this 
order.  The  solid  curves  are  the  results  of  the  tetra- 
hedron treatment  of  Sec.  IV  and  the  dot  curves  of  the 
pair  treatment  of  Sec.  II.  The  former  is  always  more 
gradual  than  the  latter  for  the  same  temperature. 

For  the  rough  measure  of  the  "thickness"  L of  the 
boundary  layer,  we  use  the  definition  used  by  Weeks  and 
Gilmer1  in  their  Eq.  (38): 

t = (l-2.v“L»l)/(l-2^.ll,1)  . (4.19) 

Note  that  the  center  of  the  boundary  is  between  v=  10  and 
11.  The  quantity  is  the  bulk  value.  This  "thick- 

ness" L is  plotted  in  Fig.  8.  Reflecting  the  fact  that  the 
tetrahedron  boundary  is  less  sharp  than  the  pair  bound- 
ary, the  L value  for  the  former  is  larger  than  the  latter. 

The  nature  of  the  phase  transitions  within  the  bound- 
ary, marked  by  crosses  in  Fig.  6,  is  one  of  the  prime 
interests  of  Weeks-Gilmer’s  paper1  and  of  the  present 
work  as  well.  In  Ref.  17  it  was  suggested  that  the  un- 
stable disordered  phase  below  T„  corresponds  to  the  AC 
boundary  and  the  stable  ordered  phase  to  the  IC  bound- 
ary. The  aAC-oIC  plotted  in  Fig.  6 indicates  that  this 
identification  almost  holds  but  is  not  exact  because  the 
tail  of  aAC-alc  extends  into  the  temperature  above  T„. 
However,  the  prediction  (C)  in  Sec.  IV  of  Ref.  17  that 
CT*c~aic  becomes  smaller  as  the  approximation  im- 
proves and  would  vanish  in  the  limit  of  the  exact  treat- 
ment is  supported  without  doubt  by  comparison  of  the 
two  aAC-alc  curves  in  Fig.  6. 

The  curves  for  "thickness  L were  studied  by  Weeks 
and  Gilmer1  to  obtain  a clue  for  the  nature  of  the  transi- 
tion point  T ,.  The  two  curves  in  Fig.  8 do  not  give  any- 
thing definite  to  answer  the  question. 

As  was  commented  in  Sec.  IV(E)  of  Ref.  17,  time- 


dependent  analysis  of  the  boundary  motion  still  remains 
a hopeful  method  of  clarifying  the  nature  of  T0. 


4.  Summary 

The  present  paper  shows  that  the  iterative  calculation 
proposed  by  Weeks  and  Gilmer1  can  be  interpreted  as  a 
modified  form  of  the  NI  calculation  of  the  cluster-varia- 
tion method  the  author  had  proposed  earlier. 5 The  rein- 
terpretation of  WG’s  method  in  the  light  of  the  NI  tech- 
nique makes  it  possible  to  extend  their  idea  and  to  cal- 
culate boundary  structures  using  larger  clusters  and 
improved  calculations.  This  is  demonstrated  in  the 
present  paper  using  the  tetrahedron  approximation  for 
the  bcc  (110)  boundary.  The  paper  also  shows  how  a 
tetrahedron  can  be  used  in  calculating  the  excess  bound- 
ary free  energy  of  the  bcc  boundary  using  the  scalar- 
product  expression. ls* " The  existence  of  the  phase 
transition  within  a phase  boundary11  is  supported  by  the 
calculation  in  the  present  paper. 
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The  scalar-product  expression  of  boundary  free  energy  for 
long-range  interacting  systemsa) 

Ryoichi  Kikuchi 

Hughes  Research  Laboratories,  Malibu,  California  90265 
(Received  5 July  1977) 

The  scalar-product  formula  of  the  excess  free  energy  cr  of  a boundary  between  two  phases  (based  on  the 
lattice  model)  is  proved  for  the  case  of  interaction  potential  of  the  range  longer  than  the  nearest  neighbor. 

The  formula  is  exp<  - A <r/kT)  = Up  'V„Vj vk)  p»'\v  u„)]>'J  exp[<*m(v,.v, v„)~ 

a""(w„w2 )],  where  A is  the  sectional  area  parallel  to  the  boundary,  o,  is  a configuration  of  an  ith 

plane  parallel  to  the  boundary,  p")(vl,...,vi)  is  the  probability  that  k consecutive  planes  in  the  bulk  1 

phase  take  configurations  v, , and  the  summation  goes  over  all  configurations  vlt...,vk  The 

variable  a",(V|,i/2>...,v»)  is  a Lagrange  multiplier  to  guarantee  continuity  of  pm: 

11  Wi.t-j •'»)=  ,"(<'i,*'i,- ■•»'*,/»)■  The  expression  is  checked  by  two  examples.  The  cr  for 

the  two-dimensional  Ising  model  is  calculated  using  a 3x2  cluster  (i.e.,  a double  square  cluster  made  of 
six  lattice  points)  with  the  "3"-side  perpendicular  to  the  boundary,  and  is  compared  with  the  previous  cr 
calculated  with  a 2x3  cluster  (with  the  “3"-side  parallel  to  the  boundary).  The  calculated  o-’s  agree  well 
when  the  a terms  are  included.  As  a second  example,  surface  tension  cr  of  a liquid  of  a two-dimensional 
lattice  gas-liquid  model  (in  which  the  first,  second,  and  third  neighbor  pairs  are  excluded,  and  the 
fourth  and  fifth  neighbor  pairs  attract)  is  calculated.  It  is  then  compared  with  cr  calculated  by  a sum 
method  (which  calculates  the  equilibrium  state  of  a sandwich  system  made  of  the  gas  and  the  liquid  phases 
with  the  boundary  between  them).  The  agreement  between  the  two  calculations  supports  the  correctness  of 
the  proposed  cr  expression. 

helped  solve8,7  the  puzzle  in  the  square  gradient  theory8 
of  the  boundary  structure,  and  was  also  useful  In  pre- 
dicting the  existence  of  a phase  transition  within  a 
boundary.  9 In  spite  of  these  successes,  problems  re- 
main with  regard  to  the  SP  formula.  We  consider  these 
problems  in  this  paper. 

When  Clayton  and  Woodbury1  proposed  the  SP  formula 
(1.1),  the  proof  was  not  as  complete  as  had  been  hoped 
for.  Further  study  by  the  author8,10  has  supplemented 
the  proof,  but  it  still  needs  improvement.  Another, 
somewhat  related,  problem  is  the  case  of  interaction  po- 
tential of  longer  range  than  the  nearest- neighbor  dis- 
tance. The  long-range  interaction  case  was  discussed 
in  Ref.  6,  Sec.  V,  but  the  reported  expression  was  in 
error. 

The  corrected  proof  of  the  SP  expression  is  presented 
in  Secs.  II  and  III  in  the  present  paper  for  the  general 
case  of  a long-range  interaction.  An  example  of  an  Ising 
model  is  presented  in  Sec.  IV  and  another  example  of  a 
gas-liquid  phase  boundary  in  Sec.  V. 

II.  PROOF  OF  THE  SP  EXPRESSION  OF  BOUNDARY 
FREE  ENERGY 

This  section  presents  the  corrected  proof  of  the  SP 
formula  of  the  boundary  free  eneigy  for  interaction  po- 
tentials of  longer  range  than  the  nearest-neighbor  dis- 
tance. In  the  proof,  we  avoid  using  the  inverse  of  the 
transfer  matrix  P,  which  was  used  in  the  proof  reported 
in  Ref.  10. 

In  the  system  being  considered,  there  is  one  phase 
boundary  between  the  left  (I)  and  the  right  (II)  phases. 

A lattice  plane  parallel  to  the  boundary  is  simply  called 
a parallel  plane.  The  position  of  a parallel  plane  is 
designated  by  i.  Each  parallel  plane  contains  N lattice 
points;  a configuration  of  the  plane  is  designated  by  a 
Greek  letter  (p,  u,  £,  or  77)  that  takes  values  from  1 


I.  INTRODUCTION 

The  purpose  of  this  paper  is  to  report  how  the  scalar- 
product  formula  of  the  boundary  free  energy  can  be  ex- 
tended to  cases  of  interactions  of  longer  range  than  the 
nearest-neighbor  distance. 

The  following  formula  of  the  boundary  excess  free  en- 
ergy at  per  unit  area  was  first  proposed  by  Clayton  and 
Woodbury1: 

(l.D 

V 

where  A is  the  sectional  area  of  the  system  parallel  to 
the  boundary,  and  0 = l/kT.  When  we  consider  a crys- 
talline plane  parallel  to  the  boundary  (called  a parallel 
plane,  for  short),  v in  (1.1)  denotes  one  of  the  configu- 
rations of  the  parallel  plane.  The  probability  that  a 
parallel  plane  takes  a configuration  v in  the  bulk  phase  I 
is  written  as  p(I>(v),  and  the  corresponding  quantity  for 
the  bulk  phase  II  is  pttl)( v).  When  we  regard  the  array 
of  p<v(v)l,z  for  v = l,  2,  ...  as  a vector,  the  expression 
(1.1)  can  be  interpreted  as  a scalar  product  of  the  two 
bulk-phase  vectors  {p<,,(l/)1/8}  and  {/><n’(i')1/8},  and  thus 
we  may  call  (1.1)  the  scalar-product  (SP)  expression  of 
the  boundary  free  energy. 

In  the  SP  formula  (1.1),  v Is  the  configuration  of  an 
infinitely  extended  parallel  plane  of  lattice  points. 
Therefore,  to  use  (1.1)  for  numerical  computations,  it 
is  necessary  to  approximate  it  using  configurations  of  a 
finite-sized  cluster.  This  approximation8, 8 was  done 
using  the  cluster-variation  method,  and  the  results  were 
checked  by  comparing  them  with  the  exact  results  of 
Onsager4  and  of  Fisher  and  Ferdinand5  for  the  two-di- 
mensional Ising  model.  The  comparison  supports  the 
correctness  of  the  SP  expression.  The  SP  formula 
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through  A*,  A'  being  the  number  of  species  in  the  sys- 
tem. For  simplicity,  the  formulation  in  this  section  is 
done  for  a system  which  does  not  need  superlattices. 

In  demonstrating  the  long-range  interaction  formula- 
tion. we  use  the  probability  function  />,(p,  i>,  £),  which 
encompasses  three  consecutive  parallel  planes  with  the 
ith  plane  at  the  center.  The  three-plane  formulation  is 
sufficient  to  Induce  the  general  case  from  it.  In  writing 


the  free  energy  F in  terms  of  />((p,  v , £),  we  use  the 
cluster-variation  (CV) approach11  and  apply  it  to  a cluster 
of  three  consecutive  infinitely  long  lines.  The  method 
is  exact  (contrary  to  a widespread  misconception  that  the 
CV  method  is  always  approximate)  since  the  system  is 
pseudo-one-dimensional.  The  method,  which  was  ex- 
plained in  Ref.  6,  Secs.  II  and  V,  leads  to  the  following 
expression  for  the  free  energy  of  the  entire  inhomoge- 
neous system: 


f Z Z «(M,  p,  £>Mm,  V,  {) 

I U , V,  t 

£ £[/>,-!/»(*!,  v)|. 2 £ £[/>,., /t(p,£)i-  £ fcfMn.i'.oi! 

* kT ^ u£  {<*i-i/,(p.  p)-  a(.w,(r.  £)}MM,  p,  £)+£ X,jl  - £ p,(M>  £)j 


(2.1) 


The  notation  differs  slightly  from  that  used  In  Ref.  6. 

In  the  first  term,  *(p,  e,  £)  is  the  energy  per  three-plane 
region;  the  £ operator  is  defined  as 

£[.r|  * x inx  - x ; (2.2) 

and  />(.,/, (p.  p)  is  the  probability  that  the  (i - 1 )th  and  the 
ith  planes  take  configurations  p and  e,  respectively. 

The  Lagrange  multiplier  a’s  are  used  to  satisfy  the  con- 
tinuity requirements: 


Pi- i/t<P.  p)  = £/>(-,<£,  P,  »’)  = £/>,( p,  p,  £)  . (2.3) 

t < 

The  sign  of  or  in  this  section  is  reversed  from  that  of 
Refs.  3 and  6.  The  last  terms  in  (2. 1)  are  for  the  nor- 
malization of  />,  for  every  i, 

£ Mm,  p,  £)  = 1 , (2.4) 

**,  ¥,( 

and  X,  are  the  Lagrange  multipliers. 

Minimizing  F in  (2. 1)  with  respect  to  />,( p,  r,  £)  yields 
Mp,  p , £)  - exp[p)X,  - 0c(p,  v,  £)| 

*[/><-i/t(P.i')/>).1/1(p,£)l,/t 

xexpfa^./.fp,  p)-  a„1/8(p,  £))  . (2.  5) 

When  F is  a minimum,  we  can  derive 


'=*•-£  £ 

< O.O.I 


8 F 

9Pi(v,  V,  £) 


=£*,  , 


(2.6) 


which  shows  that  X,  can  be  Interpreted  as  the  local  free 
energy.  For  further  transformations  It  is  convenient 
to  Introduce 


#n/i(i',  £)"[Mi/i(p,  £)]1/8  exp[-  a,.l/t(v,  £)]  , 
*t-i/i<M,  p)"(Mi/i(p,  vil’^expfan/jtp,  v)|  , (2.  7) 

T(p,  v,  £)*  exp[-  0<(p,  v , £)]  , 
and  write  (2.  5)  as 

Mp,  p,  £)  = exp p) 

*r(M,  p,  £)  . (2.8) 


I 

It  is  important  to  note  from  (2.  7)  that 

££*...  /a(M,  P)tf(.i/8(p,  p)  1,  (2.9) 

u v 

in  which  the  normalization  of  />(<1/8(p,  r)  is  used.  By 
using  the  definitions  of  f>  and  li  in  (2.7)  and  substituting 
(2.  8)  into  the  continuity  relations  in  (2.3),  the  following 
two  relations  can  be  derived: 

X'i-l/l(M,  p)  exp(0X,)£  Hp,  c,  OaVwjU'.  0 . (2. 10a) 

t 

Mt/*(P,£)  exp(0x,)£  T(p,  v,  £);,,., /t(p,  v)  . (2.10b) 

So  far  the  transformations  are  similar  to  those  in  pre- 
vious papers.8,10  Now  we  formulate  the  new  proof  of  the 
SP  expression.  We  choose  any  two  integers  i and  j with 
the  requirement 

• (2.11) 

We  first  sum  the  following  expression  over  £ and  use 
(2. 10a);  we  then  sum  over  p and  use  (2. 10b)  to  derive 

££  £ *r-i/8<M,  r)r(p,  v,  £) 

i 4 M f 

= £ £ *i-i /i(M.  p)[exp(-  0X,)x',.i/8(P,  p)| 

££  [expl-^X, )/■,., ^(r,  {)!*,., „(p,  £)  . (2.12) 

When  i=j,  use  of  (2.9)  in  the  sums  reduces  (2.12)  to 
exp(-  0\,).  In  general,  we  can  rewrite  (2. 12)  and  de- 
rive a general  recurrence  relation  for  shifting  (»’,»)  into 
(«  + l,  J + l): 

££  *1-1  /«<M,  p)tfi-i/«(p,  p) 

U V 


■ exP(~  + ^/)££*ct/l(M,  p)#r.i/t<M,  p)  . 

4 V 

(2.13) 

We  start  with  i - w and  j m + 1,  and  operate  the  re- 
currence relation  (2.13)  Af(>  2m)  times  to  obtain 


$ 


\ 


i 

i 
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EE  V ) 

U V 

= expj-  0(\^,  +X^m.l  + • • • 4 A.)+S(X„ll  + XJ„l.(  4 • . • 4 A*.,)]  5^53  v)gm.M. i/2(m,  p)  . (2. 14) 

U V 

Note  that  Amil  + • • • + A^.v-i  cancels  between  the  - 0(  ) and  the  4 £(  ) sums.  Since  M in  (2. 14)  can  be  arbitrarily 
large,  we  let  It  become  infinite.  Then  *-..*-1/2  and  gm.K.t/i  can  be  replaced  by  the  values  for  the  bulk  11  phase  so 
that 

lim  EE  v)gm* S4Wt(M,  v)  = EE  hw(n,  p)g<n>(n.  u)  = l , (2  15) 

where  we  have  used  the  normalization  (2.  9). 

We  choose  that  the  geometry  of  the  sum  in  (2. 14)  is  such  that  the  boundary  lies  between  - m and  4 m when  the  m it 
made  large.  In  the  limit  of  the  A-sum  part  leads  to 

♦m 

-01im£(A,-A  l0>)«-&4(7,  (2.16) 

where  A*0'  is  the  bulk  value  and  o is  defined  as  the  boundary  excess  free  energy  per  unit  area.  In  the  limit  as  m be- 
comes very  large,  h_m_l/z  and  gm.uz  in  (2. 14)  approach  the  hil)  and  g110  functions  for  the  respective  bulk  phases 
Therefore,  we  obtain  from  (2. 14)  for  the  limit 

^A#=£E^,W)ir,,,Vv).  (2  17) 

Writing  this  explicitly  using  (2.7),  we  obtain 

^A*=L£l/'<,’^.^<,,,(4,v)]‘>‘exp(a,“>(p,v)-a,,,>(M,v)l,  (2  i8) 

|1  V ' ' 

where  the  superscripts  I and  II  indicate  the  two  phases. 

The  proof  in  this  section  can  easily  be  generalized  to  the  case  where  the  probability  functions  p(A„  At,  . . . , A,.,) 
of  f>4l  parallel  planes  are  used  in  formulating  the  free  energy.  Then  (2.18)  is  generalized  to 

E "•  E ^,^^^.-.,^)^<,,Hv,,vft...,v,))‘»e^a',,(v,,vt,...,v,)-a‘»»(vlfvl,...,v*)l  , 

(2.19) 

I — 

in  which  the  a’s  are  the  Lagrange  multipliers  Introduced  The  proof  of  the  SP  expression  in  Sec.  V of  Ref.  6 for 
in  the  bulk  phase  to  satisfy  the  long-range  interaction  case,  which  incorrectly 

E(1)  omitted  the  a factor,  is  to  be  replaced  by  the  correct 

P d>  ■'ii  vt, . . . , vk)  proof  in  the  present  section. 


=Epm(v„  v2,  ...,v„  $),  (2.20) 

and  can  be  shown  to  obey  the  antisymmetry  relation 

„„„,) . (2.21) 

When  the  I and  II  phases  are  identical,  the  normalization 
of  p makes  a In  (2. 19)  vanish,  as  expected. 

There  is  one  remaining  puzzle  In  the  proof  of  this  sec- 
tion. We  required  the  condition  i<j  in  writing  (2. 12). 
The  choice  i = -m  and  j = m4l  in  (2.14)  satisfies  this 
condition.  However,  the  equalities  in  (2.12)  and  hence  In 
(2. 13)  seem  to  hold  even  when  i>j.  If  we  allow  i>j  in 
(2.13)  and  use  it  M times  starting  with  the  values  i=m  + l 
and  j = - m,  this  procedure  is  equivalent  to  reversing  the 
signs  in  front  of  the  A’s  In  (2. 13)  and  (2. 14).  The  sign 
reversal  results  in  the  expression  exp(4  0Ao)  on  the 
left-hand  side  of  (2. 17);  this  obviously  is  the  wrong  re- 
sult. The  interpretation  of  the  wrong  result  and  the  rea- 
son why  we  need  the  condition  i<  j in  (2. 11)  are  discussed 
in  the  Appendix. 


III.  APPLICATION  OF  THE  FINITE-SIZE  CLUSTER 
APPROXIMATION 

In  the  general  SP  expression  (2. 19),  v,  represents  one 
of  the  Ks  configurations  of  a parallel  plane  (as  was 
mentioned  in  Sec.  n),  and  the  number  K"  is  practically 
infinite  as  far  as  the  computer  calculation  is  concerned. 
Therefore,  to  make  use  of  (2. 19)  in  calculating  a.  It  Is 
necessary  to  introduce  a certain  scheme  which  reduces 
(2. 19)  into  a tractable  form.  The  case  of  the  one- plane 
probability  function  p{v)  was  done  in  Refs.  2,  3,  and  10 
using  the  finite-size  cluster  (FSC)  approximation  of  the 
cluster-variation  (CV)  technique. 

Ih  the  general  case  of  (2.19),  however,  the  reduction 
scheme  due  to  the  FSC  method  faces  a still-unsolved 
problem  of  reductng  a(,,(v,,  vt,  . . . , vk)  Into  a tractable 
form.  How  to  handle  the  a’s  is  studied  in  this  section  by 
examining  the  homogeneous  phase  version  of  (2.  5)  as  an 
example. 

To  avoid  cumbersome  notation  for  a large-size  three- 
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FIG.  1.  A 2x3  and  a 3x2  squares  used  in  the  formulation. 


dimensional  cluster,  we  present  the  discussion  in  this 
section  based  on  a two-dimensional  system.  The  paral- 
lel plane  we  considered  in  the  previous  section  now  be- 
comes a parallel  line. 


ble  II  in  Ref.  3.  This  is  not  surprising  since  the  bulk 
phase  does  not  depend  on  the  direction  in  which  the  basic 
cluster  is  placed.  The  letter  t in  Table  1(c)  stands  for 
the  ith  species  (for  example  i 1 for  a plus  spin  and  i 2 
for  a minus  spin). 

The  purpose  of  this  section  is  to  study  the  meaning 
of  a in  (3. 1).  To  do  this,  we  look  into  the  “superposi- 
tion” expressions  in  the  bulk  phase.  When  a DS  is  used 
in  the  CV  method,  the  free  energy  minimization  leads 
to  the  following  expression  [which  was  derived  in  (3. 13) 
of  Ref.  3]: 

expUo/3  + aJJJj  + + >,%  + >"„*,)  , 

(3.2a) 

where 


exp(- 


(3.2b) 


The  notations  are  slightly  different  from  those  in  Ref. 

3.  The  superscript  (I)  indicates  the  bulk  phase  (I); 

is  the  energy  per  DS;  and  exp(X|>/3)  is  the  normal- 
ization factor,  where  A0  has  the  meaning  of  the  free  en- 
ergy per  lattice  point  [in  contrast  to  (3. 1),  in  which  X(l) 
is  the  free  energy  per  N lattice  points  in  one  long  line). 


TABLE  1.  The  double-square  cluster. 


In  the  homogeneous  phase  I,  we  can  write  (2.  5)  as 

Pm(p,  v,  £)  = exp[0A(n  - fc (p,  v,  |)] 

x[/>(,,(p,  v)p(t)(v,  £)ll/2 

x exp[a(,>(p,  v)+  a(I)U,  v)] , (3.1) 

where  p,  v,  or  £ represents  one  of  K " possible  config- 
urations of  a parallel  line  made  of  N lattice  points,  and 
/>n,(p,  v,  £)  is  the  probability  that  the  three-line  cluster 
takes  the  configuration  specified  by  (p,  v,  £)  in  the  bulk 
phase  I.  We  have  used  in  writing  (3. 1)  the  antisymme- 
try relation  alv{v,  £)  = - a{,)(£,  v),  which  is  a special 
case  of  (2.  21). 

In  the  actual  computation  of  o,  we  use  an  FSC.  Our 
program  is  to  rewrite  plI)( p,  v,  £)  of  (3. 1)  In  terms  of 
quantities  of  FSC  and  examine  the  a part  of  it.  As  an 
example,  we  take  a two-dimensional  square  net  for  the 
system  and  a double  square  (DS),  illustrated  in  the  low- 
er part  of  Fig.  1,  as  the  FSC.  The  system  does  not  in- 
clude sublattices,  but  is  otherwise  general;  it  can  be 
regarded  as  an  Ising  model  system,  a disordered  phase 
of  an  alloy,  or  a lattice  model  of  a gas-liquid  system 
(the  latter  will  be  treated  in  the  next  section). 

In  Ref.  3,  we  used  the  upper  cluster  of  Fig.  1,  and 
used  the  simpler  o formula  (1. 1).  To  distinguish  the 
two  directions  of  the  clusters  In  Fig.  1,  we  will  call  the 
upper  one  a 2x3  or  a DS(n)  and  the  lower  one  a 3x2  or 
a DS(i).  Although  the  two  cases  look  similar,  the  lower 
one  needs  the  a’s,  which  are  of  concern  In  the  present 
paper. 

The  variables  we  use  in  treating  the  bulk  homoge- 
neous phase,  (I)  or  (II),  using  the  DS(l)  are  listed  in  Ta- 
ble I.  This  table  happens  to  be  exactly  the  same  as  Ta- 
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Our  aim  Is  to  write  v,  £)  ot  (3. 1)  using  h<1’  of 

(3.2)  and  then  to  find  the  quantity  corresponding  to 
a"'(p,  t')  * *»“’(£.  v)  of  (3.1). 

In  (3.2a),  a and  > are  l.agrange  multipliers,  which 
are  Introduced  to  satisfy  the  continuity  relations 

,.u>  S'  „,in  V1  „.<*»  it  <hk\ 

> <i*l  JL,  M / ■ MniHH  13.3a) 

• .« 

for  a,  and  the  rotational  symmetries 

‘ i'uii  (3.3b) 

for  >.  The  relations  In  (3.  3a)  for  a are  similar  to  the 
continuity  requirement  In  (2.3),  which  Introduces 
Of. i /flu.  r).  The  rotational  requirement  did  not  appear 
In  Sec.  If  and  needs  a special  discussion  at  the  end  of 
this  section. 

When  the  system  Is  a bulk  homogeneous  one,  ir  satis- 
fies the  symmetry  relations 

‘ K'llrtit i (3.4a) 


k’i!*!*.  ie}}|tm  . (3.  4b) 

We  can  prove  that  these  relations  (3.  4)  In  turn  lead  to 
the  following  symmetry  relations  for  or  and  y: 

Oifii  3 - . Or(%  ail,’,  , (3.5a) 


1.1H  1.1H 

>h*i  - f ft/i 


> ,%  -riH.  • 


In  Eqs.  (3. 2)—  (3.  5),  the  variables  with  the  superscript 
(1)  are  those  of  the  equilibrium  state  In  the  hulk  phase  1. 
Now  we  use  these  variables  to  write  /><l’(p.  v,  £)  In  (3. 1). 
The  variables  n,  v,  and  £ represent  any,  In  general  ex- 
cited. state  of  the  three  consecutive  lines  In  the  bulk 
phase  l.  In  designating  such  an  excited  state,  we  again 


use  the  variables  listed  in  Table  1.  but  this  time  without 
the  superscript  (1).  In  other  words,  when  we  look  at  any 
3 '2  cluster  portion  of  the  three-line  configuration 
(p,  v,  £),  the  probability  of  finding  the  configuration  i-J- 
k-1-tn-M  Is  written  as  Using  this  quantity,  we 

can  write  pal(p,  e,  £)  as 

f>(,,(p.  v,  £)  -- ]*I  [M’Jfi,.y*}!il  ••  (Ns’i/*,..) , (3.  6) 

where  a double  asterisk  Is  a FORTRAN  notation  meaning 
“raised  to  the  power  of"  and  is  used  In  this  paper  to 
avoid  subscripts  on  a superscript.  The  number  (6) 
above  the  product  sign  in  (3.  6)  Indicates  that  this  Is  a 
sixfold  product  over  i,  j,  k,  l,  w.  and  n.  In  (3.  6)  the 
factor  In  the  square  brackets  is  the  conditional  proba- 
bility of  finding  i-k-tn  next  to  j-l-n  when  the  latter  Is 
known.  The  logic  of  writing  (3.  8)  Is  the  same  as  that 
presented  In  Refs.  2,  3,  and  10. 

We  now  go  back  to  (3.2)  for  the  FSC  and  use  it  for 
(3.  6)-  The  resulting  expression  can  lie  sim- 
plified In  several  steps.  Resides  to, we  use  r’s  and 
z’s  as  shown  In  Table  1: 

m ,n 

• (3.7) 

7,l,n 

We  note  the  following  relations: 

iti 

nui&/.jB>..w*«.) 


fi(.v!iV.v5!’)**(Nw, 


Using  (3.8)  and  (3.2),  we  can  transform  /><"(p,  r.  £)  In 
(3. 6)  to 


/>'”(q.  If,  £)  exp  ( ^ dod-  PtiiHmMvilHml)  IT  ('yjf  **(  2 

* exp(V  (ojii,  • <*«ii  ♦ >i!ii  + ri'i ,)N»uumn)  . IS.  9> 

Next  we  compare  (3.9)  with  (3.1).  Rased  on  the  normalization  of  «•.  we  can  identify 

A'1’-  £ (3  >0) 

because  x"‘  Is  the  free  energy  per  .V  lattice  points  and  ,\0  Is  that  per  lattice  point.  For  the  energy,  we  can  Identify, 

«(P.  *'.<)•  . (3.11) 

Similar  to  (3.  6),  we  can  write 

(4) 

p'"(p,  v ) * U’iRi/y}*’)  •«  (JVt’i,*,) . tS.  12a) 

/>"’(v.£)«n  <3. 12b) 


(3.  12bl 


Because  ot  the  definition  ot  i>,„,  In  (3.7),  t>)(M  In  (3. 12a)  can  be  replaced  by  m , and  the  product  In  (3. 12)  can  be 
changed  over  six  indices  i,  j,  k,  I,  m , and  n as  In  (3.  9).  Thus  we  can  rewrite  (3.  9)  In  a form  close  to  that  ol  (3. 1 ): 

/><n(p,  v,  £)  exp[0\n’  - d«(p.  v.  £)1[/><,,(P.  i')/>'”(v,  £)lw* 

^ ^ • 13.13) 
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FIG.  2.  A 3 * 3 at*  u art*  used  In  the  formulation. 


When  we  compare  (3. 13)  with  (3. 1),  we  are  tempted 
to  identify 


a"\u,f)-N  £ + . (3.14) 

i 

This  identification,  however,  cannot  be  justified,  be- 
cause arjjj,  satisfies,  as  in  (3.5a),  the  same  kind  of  the 
antisymmetry  relation  as  a^Vi,  c,,  . . . , t>.)  In  (2.21), 
while  the  symmetry  property  of  yJJJ,  in  (3.  5b)  is  of  a 
different  kind. 


When  we  think  of  the  difference  between  the  formula- 
tion in  Sec.  II  based  on  the  infinitely  long  line,  and  that 
of  the  present  section  based  on  the  FSC,  we  note  that 
the  former  does  not  need  y’s  because  the  rotational  sym- 
metry is  automatically  incorporated  in  the  infinitely  long 
cluster.  This  situation  suggests  that  even  for  the  FSC 
treatment  the  effect  of  >'s  would  be  small;  this  is  exact- 
ly the  case  that  was  explicitly  pointed  out  in  Sec.  Ill  of 
Ref.  3.  In  other  words,  even  with  yllt,  0,  the  aniso- 
tropy of  the  system  is  small  and  can  be  neglected. 
Based  on  this  observation,  in  reducing  Sec.  II  to  a tract- 
able FSC  scheme  of  this  section,  we  propose  to  tolerate 
a slight  anisotropy  of  the  system  and  make 

virti  0.  (3.15) 

Thus,  instead  of  (3.14),  we  Identify 

a"\ti,v)-N  T . (3.16) 

« .77*. i 

The  result  (3. 16)  can  be  easily  generalized  to  larger 
clusters.  We  show  only  one  example.  When  we  use  a 
3x  J FSC  to  calculate  the  bulk  phase  using  the  nine-point 
variables  ic,„. based  on  Fig.  2,  am(p,  i<)  in  (3. 16) 
Is  to  be  identified  as 


aa'(n,  v) 


~U>  „<!> 

Ik,  Ilk,  Inn  ♦ 


(3.17) 


where  „„  is  the  Lagrange  multiplier  introduced  for 
the  transitional  symmetry  condition: 


„ti» 

"l/*,  Inn 


V 

tr.. 


„,ti> 

u Ilk,  fm *.<>»« 


o.P.t 


14 


,CI> 

>>Pl,Uk,  linn  • 


(3.18) 


The  FSC  expression  of  o"’(p,  v)  in  (3. 16)  or  (3. 17)  is  the 


one  which  is  to  be  used  in  calculating  exp(-  d/to)  from 

(2.181. 


IV.  TWO-DIMENSIONAL  ISING  MODEL 

As  we  did  in  our  previous  publications,  *'*  we  go  back 
to  the  two-dimensional  lsing  model  to  check  if  the  new  a 
terms  in  the  SP  expression  are  correct.  We  use  the  DS 
cluster  of  Sec.  Ill,  and  the  only  difference  betw  een  the 
present  section  and  Sec.  V of  Ref.  3 is  that  we  place  the 
DS  as  perpendicular  to  the  boundary,  as  in  the  lower 
part  of  Fig.  1.  rather  than  the  parallel  position  of  Ref. 

3 shown  in  the  upper  part  of  Fig.  1. 

The  expression  for  the  excess  boundary  free  energy 
o is  (2. 18),  combined  with  (3. 16)  for  a"’.  Using  FSC 
notation,  we  obtain 


exp(- ,Voao,nd)  £ 

t»i/*il 


"ujuM-II.v 

n(.,(vfrv!jf')..(i.vv 


l-v-W 

,) 


x <>xp[.V  ^ ^ (arJJi,  - ajJ{!)t>lrtlJ  , (4.1) 

where  o(10)  is  the  excess  free  energy  per  unit  length,  and 
a is  the  lattice  constant.  The  variables  with  super- 
scripts are  those  for  the  corresponding  bulk  phases, 
while  rutl  and  yt,  without  superscript  are  the  boundary 
variables  to  be  determined.  The  definitions  are  in  Ta- 
ble I.  The  D factor  is  the  weight  factor  given  as 


for  which  we  have  the  relations 


(4.2) 


(4.  3' 


For  the  derivations  of  (4. 1)  and  (4.  2),  Refs.  2 and  3 may 
be  referred  to.  As  we  have  been  discussing  so  far,  the 
new  aspect  in  this  paper  is  the  last  factor  In  (4. 1)  con- 
taining a. 


To  calculate  a(10)  from  (4. 1),  it  Is  convenient  to  find 
the  maximum  of  the  logarithm  of  the  summand  on  the 
right-hand  side,  and  in  so  doing  use  the  natural  itera- 
tion (Nl)  method. 5,10  The  Iteration  converges  with  ease 
and  the  results  are  shown  in  Fig.  3 by  the  solid  curve 
marked  double  square  (1).  It  is  satisfying  to  see  that 
DS(l)  is  close  to  DS(n),  which  is  the  previous  result 
reported  in  Fig.  3 of  Ref.  3.  It  is  particularly  note- 
worthy that,  when  we  deliberately  put  oj'i,  aJJJJ  0 in 
(4. 1),  the  a curve  moves  way  up  to  the  dotted  position 
In  Fig.  3;  this  fact  shows  the  correctness  of  the  a terms 
in  (4. 1)  and  hence  in  the  general  formula  (2. 19). 

Because  of  the  antisymmetry  property  of  aj'J,  In 
(3.  5a),  we  ran  replace  the  last  exponential  factor  con- 
taining a in  (4. 1 ) by 

ro8h[-v  ( - «;!{;)(>„„]  . 14.  4' 

Therefore  the  a curve  for  which  a is  taken  into  account 
is  always  lower  than  the  one  for  which  a is  left  out.  The 
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kT  It 

FIG.  3.  Excess  free  energy  a at  the  boundary  of  square  net 
Ising  model. 


DS(1)  curve  and  the  dotted  curve  in  Fig.  3 satisfy  this 
general  property. 

V.  LATTICE  GAS-LIQUID  SYSTEM 

Although  the  Ising  model  example  in  the  previous  sec- 
tion is  useful  in  supporting  the  validity  of  the  a factor 
in  (2. 18)  and  of  the  expression  (3. 16),  the  quantity  o for 
the  Ising  model  as  shown  in  Fig.  3 can  be  calculated  by 
always  starting  with  a double- line  cluster  parallel  to  the 
boundary  and  applying  the  formula  (1.1).  In  other 
words,  when  the  interaction  is  of  the  nearest- neighbor 
type,  we  can  do  without  the  new  formula  (2. 18)  or 
(2. 19).  The  new  formula  is  needed  only  when  the  inter- 
action becomes  longer  range.  We  show  an  example  of 
the  latter  case  in  this  section. 

We  calculate  the  surface  tension  of  a two-dimensional 
gas-liquid  phase  boundary  using  a model  proposed  by 
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FIG.  4.  Interacting  pairs  AB  and  AC. 


ulf 

FIG.  5.  Grand  potential  G versus  chemical  potential  near  the 
coexistence  point  for  kT/ f =0.50. 


Orban  et  al.li  The  model  is  illustrated  in  Fig.  4.  When 
an  atom  sits  at  A,  no  atoms  can  come  to  the  first,  sec- 
ond, or  third  neighbor  points  (indicated  by  crosses). 

An  atomic  pair  AB  at  the  fourth  neighbor  distance  at- 
tracts with  the  potential  - gt  («>  0),  and  a pair  AC  at  the 
fifth  neighbor  attracts  with  - t.  In  their  example,  Or- 
ban et  al.  chose  g = 1.2,  and  we  will  use  the  same  num- 
ber. They  were  interested  in  deriving  three  phases, 
gas,  liquid,  and  solid,  but  in  the  present  paper  we  will 
discuss  the  gas-liquid  transition  only.  We  will  use  a 
3x3  cluster  first  to  calculate  homogeneous  liquid  and 
gas  phases,  and  use  the  results  to  write  (2. 18)  to  eval- 
uate o. 

Since  the  formulation  of  the  bulk  phase  calculations 
follows  the  standard  CV  method,  only  some  important 
points  are  discussed  in  this  section.  In  working  with 
gas  and  liquid  phases,  we  fix  the  temperature  and  the 
chemical  potential  p.  The  advantage  of  fixing  p rather 
than  composition  in  the  treatment  of  phase  diagrams  has 
been  previously  discussed. 15-15  When  p is  fixed,  the 
thermodynamic  potential  which  is  minimized  is  not  the 
free  energy  F = E - TS  but  the  quantity 

G = F-pNA  (5.1) 

which  we  call  the  grand  potential;  JV„  is  the  number  of 
atoms  in  a system. 

Keeping  T and  p fixed,  we  minimize  G with  respect  to 
the  six  cluster  variables  for  a 3x3  cluster  which  speci- 
fy the  state  of  the  system,  using  the  natural  interation 
method.1*  The  resulting  G,  which  is  now  a function  of 
T and  p,  is  made  of  two  branches,  one  for  the  gas  phase 
and  the  other  for  the  liquid  phase  as  shown  in  Fig.  5. 

The  point  at  which  the  two  branches  cross  gives  the 
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FIG.  8.  Phase  separation  diagram  of  the  gas  and  liquid  phases. 


5 also  gives  information  about  the  densities  of  the  coex- 
isting gas  and  liquid  phases.  The  curve  marked  3^3  in 
Fig.  8 shows  the  result. 


for  which  two  phases  coexist.  The  curve  marked  by  3 x 3 
in  Fig.  6 shows  the  Tvs  Me-ti  relations. 

Since  the  chemical  potential  p is  the  Gibbs  potential 
per  atom,  thermodynamics  tells  us  that  G defined  in 
(5. 1 ) is  equal  to 

G - - f>A  , (5.2) 

where  f>  is  the  pressure  and  A is  the  area  of  the  system. 
(In  a three-dimensional  system,  the  corresponding  equa- 
tion is  G_  - f>V,  V being  the  volume  of  the  system. ) The 
value  of  G corresponding  to  pc  gives  the  pressure  of 
coexistence,  The  is  plotted  against  T as  the 

3y  3 curve  in  Fig.  7. 

The  crossing  point  (i.e. , the  coexisting  point)  in  Fig. 


0.010 


0.1 


FIG.  7.  T vs  pressure  at  coexistence  for  the  gas-liquid 
model. 


A comparison  of  Figs.  6,  7.  and  8 with  the  results 
reported  In  Sec.  IV  of  Ref.  12  shows  them  to  be  in  good 
agreement  except  for  the  estimate  of  the  critical  point. 
Our  Fig.  8 shows  clearly  that  the  maximum  of  the  curve 
is  at  kT/(  0.  54,  while  Fig.  9 of  Ref.  12  estimates  that 
the  critical  value  of  kT/t  is  somewhere  beyond  0.  7. 
Later  we  discuss  an  additional  evidence  which  shows 
that  the  estimate  of  Ref.  12  is  too  large. 


¥ ai 


FIG.  9,  Surface  tension  (T  of  the  gas— liquid  system  calculated 
in  Sec.  V. 
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After  the  coexisting  phases  are  thus  determined,  we 
use  (2. 18)  to  calculate  the  surface  tension  a.  Since  (as 
we  said  at  the  outset  of  this  section)  the  thermodynamic 
potential  is  G rather  than  F,  the  quantity  a is  not  the  ex- 
cess free  energy  but  rather  is  the  excess  grand  poten- 
tial.1T,la  With  this  understanding,  we  call  o the  surface 
tension.  The  result  of  the  calculation  is  plbtted  by  the 
solid  curve  marked  SP,  3x3  in  Fig.  9.  The  dotted 
curve  accompanying  it  is  the  one  for  which  a in  (2. 16) 
is  deliberately  put  equal  to  aero.  Different  from  the 
Fig.  3 case,  the  solid  and  dotted  curves  do  not  differ 
much,  although  the  dotted  curve  (a  =0)  is  higher  than 
the  solid  curve  (a  *0),  in  agreement  with  the  general 
requirement  mentioned  at  the  end  of  Sec.  IV. 

The  reason  why  the  effect  of  a is  large  in  the  Ising 
model  (Fig.  3)  and  is  practically  nil  in  the  gas-liquid 
surface  tension  (Fig.  9)  can  be  traced  in  mathematics, 
but  the  physical  reason  is  still  to  be  determined. 

For  the  two-dimensional  Ising  model,  the  result  of  o 
calculated  by  the  SP  method  can  be  compared  with  the 
exact  calculation  due  to  Onsager,  * as  we  did  in  Fig.  3. 
Since  in  the  gas- liquid  model  of  this  section  there  is  no 
exact  calculation  to  compare  with,  we  calculate  a using 
two  different  methods  to  check  the  accuracy  of  the  SP 
method.  One  of  them  is  the  "sum"  method.  We  mini- 
mize the  grand  potential  G of  an  inhomogeneous  system 
that  includes  the  gas-liquid  boundary.  We  apply  the  CV 
method  using  a 3x3 square  of  Fig.  2 as  the  basic  cluster. 
The  technique  is  similar  to  the  one  used  by  Cahn  and 
Kikuchi17,1*  and  modified  later, 10  taking  into  account  the 
Weeks  and  Gilmer  technique. 19  In  calculation,  we  used 
80  lattice  lines  and  imposed  the  conditions  that  the  left 
two  end  lines  are  in  the  gas  phase  and  the  right  two  end 
lines  are  in  the  liquid  phase.  In  minimizing  G,  the  com- 
bination of  T and  p are  fixed,  the  latter  being  the  value 
pc->_  for  which  the  gas  and  liquid  phases  coexist  at  that 
temperature,  as  determined  in  Figs.  5 and  6. 

The  results  of  <j  calculated  by  the  sum  method  are 
plotted  by  the  broken  curve  in  Fig.  9.  When  we  com- 


FIG.  10.  Density  profile  versus  the  boundary  calculated  by  the 
sum  method. 


lattice  line  across  the  boundary 

FIG.  11.  Local  excess  grand  potential  across  the  boundary 
calculated  by  the  sum  method. 


pare  the  SP  (3x3)  curve  and  the  sum  curve,  the  former 
is  higher  than  the  latter  for  the  region  of  physical  sig- 
nificance, kT/(>  0.  4.  (For  the  temperatures  below 
about  0. 4,  the  solid  phase  is  more  stable,  as  was  shown 
in  Ref.  12,  and  the  curves  in  Fig.  9 lose  their  meaning.) 
This  situation  qualitatively  agrees  with  the  bcc  ( 110) 
boundary  results  reported  in  Fig.  6 of  Ref.  10;  in  the 
bcc  case  also,  the  SP  curve  is  higher  than  correspond- 
ing sum  curve. 

We  can  then  almost  say  that  the  SP  expression  (2. 18) 
is  supported  by  these  calculations.  There  is,  however, 
one  bothersome  feature  of  the  SP  (3x  3)  curves  in  Fig. 

9:  they  bend  over  around  kT/e  = 0. 25.  We  can  disre- 
gard, if  we  want,  this  bendlng-over  behavior,  because 
it  occurs  in  the  (T,  p)  region  in  which  the  gas  and  liquid 
phases  are  not  stable,  as  we  discussed  above.  (In  this 
regard,  we  may  quote  Barker20  who  pointed  out  that  the 
CV  calculation  can  sometimes  give  unphysical  results 
in  the  region  where  the  phase  being  calculated  is  not  the 
most  stable  one. ) However,  to  further  verify  the  SP 
method,  we  calculated  cr  using  a 3x4  cluster.  The  re- 
sults are  plotted  by  the  curves  marked  3x4,  SP  in  Figs. 
6-9.  The  solid  curve  in  Fig.  9 is  the  or  * 0 case  and  the 
dotted  one  is  the  a-0  case;  the  latter  is  higher,  In 
agreement  with  the  requirement  at  the  end  of  Sec.  IV. 
The  3x4,  SP  curve  does  not  show  the  bending  over  that 
the  3x3  curve  does,  and  behaves  more  normally;  there- 
fore, we  do  not  need  to  be  concerned  about  the  somewhat 
bothersome  shape  of  the  3x3  curve  for  low  tempera- 
tures. 

One  other  result  of  the  3x4  cluster  calculation  worth 
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pointing  out  is  the  liquid-gas  diagram  in  Fig.  8.  The 
diagram  shows  that  the  critical  temperature  calculated 
using  the  3*4  cluster  Is  bTU  0.  51.  which  is  lower  than 
the  value  lor  the  3*3  cluster.  It  is  generally  accepted 
that  the  critical  temperature  decreases  as  the  approxi- 
mation improves,  and  Fig.  8 agrees  with  this  general 
property.  This  property  is  the  other  evidence  that  the 
estimate  o 1 the  critical  region  in  Rel.  12  is  too  high. 
Although  the  surlace  tension  <j  in  Fig.  9 is  for  a two- 
dimensional  system,  the  general  feature  agrees  with 
that  of  real  three-dimensional  liquid. 

Although  they  are  not  the  prime  interest  of  this  paper, 
Figs.  10  and  11  show  the  boundary  profile  and  the  local 
excess  grand  potential  derived  from  the  sum  method. 

The  oa/t  plotted  by  the  broken  curve  in  Fig.  9 is  the 
sum  of  the  local  values  in  Fig.  11.  A part  of  the  curve 
in  Fig.  11  becomes  negative.  This  negative  part  is  con- 
trary to  what  we  expect  in  the  "central  hump"  reasoning 
used  in  the  square  gradient  treatment*  of  a,  as  was  dis- 
cussed in  Ref.  6.  and  is  worth  further  attention. 

1 


APPENDIX:  SUPPLEMENT  TO  THE  PROOF  IN  SEC.  II 

In  this  Appendix  we  first  show  the  meaning  of  the  expression  v)#(.1/t(p,  v)  which  appears  in  (2. 14),  with 

a particular  emphasis  on  why  i<  j in  (2. 11)  is  required.  Since  the  interaction  potential  in  this  formulation  is  defined 
to  be  extended  to  the  second  neighboring  plane  but  not  farther,  the  joint  probability  distribution  of  four  plane  con- 
figuration p,. i/j(m.  v,  £,  tj)  centered  at  the  position  i*l/2  is  written  without  approximation  as 

Pi. V,  £,  p)  f>iU.  V,  OPeifv,  £,  p)/pi.i/i(v,  £)  . lAl) 

We  use  (2.  8)  for  p,(n,  v,  £)  and  divide  both  sides  of  (Al)  by  r(p,  i\  £)r(v,  £,  ij)  to  obtain 

exp(£A|  +£A1.i)A(.|/i(m,  e)g(.i/|(£.tj)  = f>i.i/j(P.  v,  £,  rj)exp(d<(u,  v.  0 * $t{v,  £,  t))1  . (A2) 

We  can  continue  multiplying  the  conditional  probabilities,  as  was  done  in  (Al),  m times  to  arrive  at  the  probability 
distribution  of  (m  *l)-plane  configurations  iq,  iq,  . . . , e,  on  planes  i,  «♦  1,  ....  m as 

explSX,.,  + 0X,.,+  • ••  +3V»-i)Vi/t(i'o.  l'iVgt.«-i/i(v-_1,  v„) 

= Pii Oexpflfcleo,  e,,  !'«)♦•  •*  ♦ ^*<*'.-1.  >'■-».  ‘'Jl  • lA3) 

This  is  the  meaning  of  the  factor  h,.ut (p,  i')g,.,/,(p.  v).  In  order  that  the  expression  iA3)  be  meaningful,  m must  be 
larger  than  or  equal  to  unity;  this  is  the  condition  exactly  equivalent  to  the  condition  i<  j required  in  (2.  11)  to  accom- 
pany (2. 13). 

Since  the  left-hand  side  of  (A3)  does  not  contain  iq,  the  right-hand  side  is  also  independent  of  these 

configurations  although  they  appear  in  the  expression.  Now  we  start  with  the  expression 

Ssex p(0X(.,  + ---  + 0 v)r(p,  v,  Of!im.i,t^v.  £)  , (A4) 

* * i 

and  sum  this  first  over  £ and  later  over  p.  With  the  use  of  (2. 10).  we  obtain 
S = exp(0X,„  ♦•  ••  + 0X, ,*.,)£  , v)gi^-i/j(ix.  v) 

U V 

= exp(0X,., + • ••  +0X,.m)]jn£j!(.j/t(i',  £)gf.m.i/»(i',  £)  • (A5) 

•>  < 

The  equality  in  (A5)  shows  that  this  sum  is  independent  of  the  location  i.  When  we  choose  i such  that  the  l egion 
from  i to  i * m covers  the  boundary  (this  assignment  is  always  possible  because  the  boundary  profile  has  been  solved 
beforehand,  and  m is  to  be  made  very  large),  the  invariance  in  (A5)  leads  to  the  o expression  (2. 17). 

Note  that  in  this  proof  the  expression  S in  (A5)  Is  exactly  the  sum  of  (A3)  over  p = e0  i/„.,  and  v = iq  r„.  and  hence 

that  the  boundary  free  energy  is  closely  tied  to  the  multiplane  distribution  function  fq i i-m/fo'  Vj, . . . , i»„)  across 

the  boundary. 

Now  we  go  to  the  puzzle  mentioned  at  the  end  of  Sec.  II.  The  key  to  solve  the  puzzle  is  the  fact  that  the  probabil- 
ity expression  for  m + \ planes  ptl  iq*  e in  (A3)  requires  that  the  h function  lies  always  on  the  left  of 


VI.  SUMMARY 

The  scalar-product  (SP)  expression  of  the  excess 
boundary  free  energy  <r  writes  a in  terms  of  the  proper- 
ties of  two  bulk  phases  that  meet  at  the  boundary.  Anoth- 
er method  of  calculating  o.  called  the  sum  method  in 
this  paper,  evaluates  o as  the  difference  between  the 
free  energy  of  an  inhomogeneous  system  including  a 
boundary  and  that  of  the  homogeneous  phase.  Compared 
with  the  sum  method,  the  SP  method  has  the  advantage 
that  it  can  calculate  o more  easily  and  can  derive  some 
of  the  properties  of  o more  accurately.  In  the  papers 
published  so  far.  the  SP  method  has  been  applied  only  to 
cases  of  nearest-neighbor  interaction.  The  present 
paper  extends  the  SP  method  to  cases  of  longer  range 
interactions  than  the  nearest  neighbor.  The  extended 
SP  formula  is  in  (2. 18)  and  more  generally  in  (2. 19). 

As  examples,  o’s  are  calculated  for  a two-dimensional 
Ising  model  (nearest-neighbor  interaction)  and  for  a two- 
dimensional  gas-liquid  model  taking  into  account  up  to 
the  fifth  neighbors. 
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the  g function.  As  a mathematical  expression,  however,  we  ran  interchange  location  i<l  2 and  1 2 for  li 

and  g and  can  work  with  K|.i/illV  0-  We  now  examine  the  physical  meaning  of  this  expression.  Our 

interpretation  vs  that  this  new  formulation  corresponds  to  a different  expression  />[( im,j(Foi  vu  • • • , P„)  for  a mod- 

el of  different  interaction  energies.  We  use  a caret  to  indicate  a function  in  the  redefined  system.  The  li  and  g 
functions  are  defined  using  the  quantities  in  the  old  system  as 

v)  [pi-1/t^.  exp[-  p)|  , 

4? 0 * A<«f /a(v»  £ ) Ol1  /2expl  * a*  «i  /j(p»  0]  • (A6) 

The  energy  factor  TIm,  v,  £)  in  (2.  7)  are  replaced  by  T which  obeys 


T.  f (p,  t>,  {)r(p,  i>,  n)  = 0,  . , (A7a) 

f (p,  »»,  £)  - f(£,  v,  u)  . (A7b) 

Then  f>t(u.  v,  £)  in  (2.  5)  is  replaced  by 

e,£)  exp(-  v)t (p,  v,  Vgi.i/tW,  £)  . (A8) 

Note  that  3\,  in  (2.5)  is  replaced  by  - 3X(  in  the  new  expression  (A8).  The  (>n  ♦ 1 (-plane  probability  corresponding 
to  (A3)  is  written  as 

hi i-iU'o  ‘'i * Jit  too.  V\,  v2,  !•,)••  • f(v„.2,  c„.„  Vm)]-1 

exp(-  3^1.1“  - 3*1 -»-l  vOgl^n-l/ltom-l,  vm ) 

exp(- 3X(.!  — • • • - pX(^,.l)^,.l^2(e0,  i'l)/i,^,.1/2(e„.1,  e„)  . (A9) 

Note  that  the  h function  appears  on  the  right  of  g,  although  li  appears  on  the  left  of  g as  required.  From  the  in- 
variance argument  similar  to  (A5),  we  can  prove  that  the  boundary  free  energy  a for  this  new  system  is 


Ad  - ^U,-X'“’)<0  (A10) 

i 

as  we  said  at  the  end  of  Sec.  II.  The  negative  8 means  that  the  newly  defined  system  does  not  sustain  a stable  phase 
boundary  in  it. 


J 


.Vote  added  in  proof:  As  the  answer  to  the  comment 
at  the  end  of  Sec.  V,  Professor  J.  W.  Cahn  told  the 
author  that  the  negative  part  of  the  local  excess  grand 
potential  in  Fig.  11  can  be  understood  from  the  square 
gradient  theory  by  a partial  integration. 
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ABSTRACT 


The  ground  state  of  a binary  fee  alloy  with  nearest-neighbor 
interactions  and  tetrahedral  multiatom  interactions  (characterized  by 
two  parameters  a and  p)  is  derived  and  presented  in  a chart  which  shows 
classification  of  phases  in  the  space  of  a and  p.  Where  non-stoichio- 
metric  phases  can  occur,  the  phase  boundaries  at  T=0  are  computed  using 
the  tetrahedron  approximation  of  the  cluster  variation  method.  Duality 
of  the  composition  and  the  chemical  potential  is  emphasized,  and  results 
are  presented  and/or  interpreted  in  the  dual  approaches. 


1. 


Introduct  Lou 


Statistical  recli.mic.il  calculations  of  alloy  problems  usua  iiy 
assume  a model  with  pair-wise  interaction  energies.  Even  the  simplest 
model  with  near  neighbor  pair-wise  interactions  on  an  undistorted 
lattice  poses  a formidable  problem  which  has  not  been  solved  in  three 
dimensions.^  The  cluster  variation  method  (CVM)  has  been  a widely  used 
tool  for  performing  approximate  but  increasingly  accurate  calculations 
on  such  models. ^ 

To  make  the  models  more  realistic,  several  routes  are  ope’  . Higher- 

3 

neighbor  interactions  can  be  considered,  but  in  the  CVM  the  .plexity 
of  the  calculation  increases  with  cluster  size  and  the  minimum  cluster 
size  is  dictated  by  the  largest  interaction  distance.  An  alternate 
possibility  is  to  assume  multiatcn  interactions  such  that  the  energies 
of  clusters  are  given  by  numbers  that  cannot  be  obtained  by  summing 
pair-wise  energies.* 

These  possibilities  have  been  explored  in  fee  ordering  reactions. 

With  near  neighbor  pair-wise  interactions  the  Bragg-Williams  approximation 

gives  a completely  unrealistic  phase  diagram  which  is  unaltered  by 

4 

considering  higher  neighbor  pair-vise  energy.  In  the  CVM  the  pair 
approximation  is  unrealistic  for  different  reasons. It  is  not  until 
the  tctrnhe iron  approximation  that  a phase  diagram  which  resembles  a 
symmetric  version  of  the  Cu-Au  diagram  is  obtained.  With  practically  no 


* It  may  re  remarked  that  composition-dependent  pair-wise  energies  is 

i 

a multi  a;:- a interaction  concept  ia  which  the  energy  of  a pair  depends  on 
all  th*.  ms  in  the  volume  element  over  which  composition  is  measured. 


► 


f 

i 


I 


- A 


increase  in  computational  complexity,  fotir-.itc.u  intctaction  pa;  winters 
can  be  introduced  in  the  CVH  calculations  that  match  the  asymmetrii  Cu- 
Au  phase  diagram. ^ 

There  are  ways  of  demonstrating  that  pair-wise  interactions  are 

g 

inadequate  to  describe  alloys.  Furthermore,  there  are  quantum  mechanical 

methods  that  lend  themselves  to  the  direct  calculation  of  mult iatomi r 
9 

cluster  energies,  a that  eventually  these  may  become  available  as 

input  for  the  statistical  mechanical  calculation. 

Because  of  the  use  of  multiatom  forces  in  statistical  mechanical 

phase  diagram  calculations  in  fee,  we  undertook  a calculation^  of 

antiphase  aud  interphase  boundaries  (APB  and  1PB)  in  such  a system.  As 

11-13 

in  previous  calculations,  it  became  apparent  that  the  ground 

14  15 

state  ’ was  a clue  to  some  of  the  low  temperature  behavior.  We 
therefore  undertook  a study  of  the  ground  state  reported  in  this  paper. 

The  IPB  and  APB  arc  the  subject  of  a companion  paper. ^ 

2.  Ground  State  Energy 

We  use  the  linear  programming  method  of  Allen  and  Cahn. ^ The 
problem  is  to  minimize  the  energy  used  by  Kikuchi  and  deFontaine^ 


F./N  = 3w(l  + «)  Z2  + 4w7.2  + 3w  (1  + p)  Z ^ (1) 


where  N the  total  number  of  fee  lattice  sites,  Z is  the  fraction 

’ n 

of  fc_  - - tem  tetrahedra  containing  exactly  n B atoms  and  (4-n)  A 


atoms,!  w is  an  interaction  energy  and  u and  p are  dimensionless  uumbe 
that  express  the  strength  of  the  £:ur-body  forces.  If  pair-vise  near- 
neighbor  interactions  suffice  to  describe  the  energy,  e = p = 0. 
Equation  (1)  is  subject  to  two  constraints: 

1 = ZQ  + Zx  * Z2  + Z3  + Z4  (2) 

and  the  composition,  given  as  a fraction  x of  atoms  that  are  B, 

4x  = Zl  + 2Z2  + 3Z3  + 4Z4  (3) 

Using  the  constraints  (2)  and  (3)  to  eliminate  two  of  the  five  Z 's  we 

n 

obtain  a linear  equation  for  E in  terms  of  the  remaining  three  which 
then  is  subject  only  to  the  constraints  that  0 < Z < 1.  If  the 
coefficients  of  these  remaining  Z's  in  these  expressions  are 
positive,  E is  a minimum  when  these  Z's  are  zero.  By  successively 


t In  the  notation  of  Reference  7 and  in  a later  section  of  this 
paper  various  tetrahedral  clusters  that  differ  only  by  rotation 
are  distinguished,  and  A's  are  called  1 and  B's  are  called  2. 
Thus,  z. . , „ and  z1011  both  represent  the  tetrahedron  A0B  with 
the  B atom  on  different  corners  of  the  tetrahedron.  The  Z's 
in  the  pretent  section  are  given  in  terms  of  z's  with  four 
subscripts  by  the  following  relations 

Z0  = Z1 1 11 

“1  = *2111  *1211  + Z 1 1 2 1 + Z1 112  * etC’ 

Also,  , in  Reference  7 eq.  (4.2)  is  written  as  w in  the  present 


paper. 


eliminating  all  .nrs  of  Z's  the  minima  in  F.  are  explored.  for  i astatine, 
cq.  (1)  is  already  in  the  form  of  ZQ  and  eliminated.  Therefore,  if 

w>0,  l+u>0,  l+p>0,  (41 

then  the  zuaimi_ra  in  £ occurs  at 

E = 0 

Z1  = Z2  = Z3  = °* 

Zq  = 1 - <**,  (5) 

Z4  = *x. 

From  the  constraint  that  the  Z's  lie  between  0 and  1 it  follows  that 
0 < x < 1/4. 

If  instead  we  solve  eqs.  (2)  and  (3)  for  Z ^ and  2 

Zl  = (2-4x)  - 2ZQ  + Z3  + 2Za 
Z2  = (4x-l)  + ZQ  - 2Z3  - 3Z4 

and  use  these  to  eliminate  Z^  and  Z2  from  eq.  (1),  we  obtain 

E/N  - 3w(l+a)  (2-4x)  - 4w(4x-l) 

= - Iw^i+aa)  ZQ  - w(2-3(a+p))  Z3  - 6w(l-a)  Z,  (6) 


Thus,  when 


and 


w(l+3or)  <0 
w(2-3(«+p))  <0 


(7) 


t 

w(l-«)  <0 


t lie  nr.,  a-  vain-  foi  F.  is  given  by 


F./.S'  = 3w(l+a)  (2-Ax)  + 4w(4x-l) 


ZQ  = 7-3  = Z4  = 0,  Zx  = 2-4x,  Z2  = 4x-l,  0.25  < x < 0.50. 


These  equations  and  inequalities  are  more  easily  obtained  by 
considering  eqs . (1-3)  as  three  simultaneous  equations  in  six  unknowns. 

..  . , , 16,17 

Lsrng  Cramer  s rule,  we  obtain 


| Enm|  = ^ | knm|  Z^ 


where  j Ears | and  |knm|  are  determinants  formed  by  three  columns  from 
the  matrix 


0 3w(l+ci)  4w  3w(l+p)  0 E/N 

11  1111 
0 1 2 3 4 4x 


labelled  E,n,m,  and  k,n,m  respectively.  The  first  five  columns  of 
the  matrix,  labelled  0 to  4 are  the  coefficients  of  Z^  to  Z^  respectively 
in  eqs.  (1-3),  and  the  last  column  labelled  E is  composed  of  the  remain- 
ing terms.  Because  a determinant  in  which  a column  is  repeated  is  zero, 
the  tsrr.i  involving  and  Z^  vanish  in  eq.  (9a).  Setting  |Enm|  to 
zero  giver,  the  ground  state  energy  when  tetrahedra  n and  m are  present, 
while  tr.2  inequalities  are  given  by  requirements  that  |knm|  have  the  same 
sigr.  as  'r.-n).  Thus,  eq.  (4)  may  b^  written  in  terms  of  determinants 
J 104  , , | 3C4J  , (5)  in  terms  of  | F.04|  , (7)  in  terms  of  j 0 1 2 1 , 

j 3 12 , , and  (8)  in  terms  of  |F.12|,  where  each  symbol  represents  a 

column  n (*b). 


The  results  are  summarized  in  Tables  1 and  2,  and  in  figure  1.  Tht 
conditions  in  Table  1 are  constructed  from  the  above  rules,  bearing  in 
mind  that  permutation  of  two  columns  changes  the  sign  of  a determinant. 
Each  line  segment  in  figure  1 is  a boundary  where  one  of  the  determinants 
is  zero  and  thus  changes  sign.  Regions  iu  which  a particular  intermediate 
phase  will  appear  at  the  appropriate  composition  form  polygons  in  figure 


3.  Ground  State  Degeneracy 

This  linear  programming  method  gives  the  combination  of  clusters 
that  would  give  the  lowest  energy.  At  most,  two  types  of  clusters  can 
be  in  the  ground  state.  The  presence  of  any  others  would  raise  the 
energy  above  the  minimum.  For  example,  eq.  (6)  indicates  that  and 
can  be  in  the  ground  state,  but  the  presence  of  any  other  clusters 
raise  the  energy  above  the  minimum.  Since  each  cluster  type  by  itself 
gives  a stoichiometric  phase,  a two-phase  mixture  of  stoichiometric 
phases  would  have  the  ground  state  energy  apart  from  any  excess  due  to 
unwanted  clusters  at  the  interfaces  between  the  stoichiometric  phases. 

In  this  model,  however,  the  ground  state  can  have  many  different 
configurations.  There  are  three  sources  of  degeneracy  in  the  ground 
state.  Each  stoichiometric  ordered  phase  can  have  one  dimensional 
disorder  without  raising  the  energy.  Parallel  planes  of  certain  anti- 
phase bcundaries  (APB)  can  be  created  without  changing  the  cluster  type. 
Both  the  Ll^  (Cu^Au)  structure  and  the  DO^  (Ni^V)  structure  are  in:  de  up 
of  only  A.B  clusters.  The  D(>22  structure  can  be  thought  of  as  with 
evenly  spaced  (001)  APB's,  and  viceiversa.  In  fact,  any  distribution  of 
parallel  (201)  APB  has  the  ground  state  energy.  The  same  holds  true  for 
APB's  in  LI,  (CuAu)  of  the  type  that  converts  it  into  the  CuAu  II  struc- 


n 1 1 I 1 

FIG.  4.  Interacting  pairs  AB  and  AC, 


-j  -*•**»■  ••  — r 

method.1*  The  resulting  G,  which  is  now  a function  of 
T and  fi,  is  made  of  two  branches,  one  for  the  gas  phase 
and  the  other  lor  the  liquid  phase  as  shown  in  Fig.  5. 

The  point  at  which  the  two  branches  cross  gives  the  p 
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When  a system  is  exactly  at  a stoichiometric  composition  where  only 
one  cluster  is  present,  the  APB's  discussed  above  are  the  only  degeneracy . 
APB's  with  other  orientations  create  clusters  that  raise  the  energy.  In 
non-stoichiometric  systems,  when  the  pairs  of  clusters  present,  are  um 
adjacent  in  composition,  |ra-n|/l,  the  ground  state  phases  cannot  deviate 
from  stoichiometry  for  that  would  always  create  clusters  adjacent  in 
composition  and  raise  the  energy.  All  interphase  boundaries  (IPB)  also 
raise  the  energy,  as  do  APB's  except  the  type  discussed  above.  Thus, 
for  noil-adjacent  clusters  the  only  degeneracy  is  that  due  to  one  kind  of 
parallel  APB's. 

For  adjacent  clusters,  | m~n| = 1,  in  a non-stoichiometric  system  there 
are  several  more  sources  of  degeneracy.  Adjacent  clusters  can  be  mixed 
to  give  non-stoichiometric  phases,  nonuniform  phases  including  mixtures 
of  region  with  differing  composition,  "IPB’s"  where  such  regions  meet, 
and  APB's.  Any  such  arrangement  will  have  the  ground  state  energy  as 
long  as  only  the  two  types  of  clusters  are  used.  "Two-phase"  disper- 
sions on  any  scale  down  to  the  atomic  is  permitted. 

As  the  composition  varies  from  one  stoichiometry  to  the  next  a 
symmetry  change  is  occurring  and  there  must  be  a phase  transition. 

Thus,  as  B is  added  to  pure  A in  a (0-1)  case  the  ground  state  requires 
only  that  r.o  two  B’s  are  neighbors.  When  x is  small,  this  must  look 
like  a dilute  almost-random  solid  solution,  but  as  x approaches  0.25,  it 
must  tend  toward  an  ordered  alloy.  The  question  of  the  nature  of  the 
ordering  transition  is  examined  in  the  next  section. 

The  APB's  involved  in  one-dimensional  disorder  have  zero  excess 
energy,  f..  do  APB's  and  IPB's  when  adjacent  clusters  are  present  in  the 
ground  state  (c.g.,  cluster  pairs  0-1  and  1-2).  For  all  the  other 

there  boundaries  have  finite  energies  at  T=0.  For  the  APB's 


conditions 


KT/« 


»™el7‘  T WS  pieSSure  at  ooexiaWnoe  for  tht*  gaa-liquid 


aj 

KT/« 


o.e 


in  Sec.  V. 
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FIG.  !».  Surface  tension  cr  of  the  gas-liquid  system  calculated 


this  situation  is  similar  to  the  results  of  an  investigation  of  the  Csi  '. 
structure,*1  and  will  be  used  to  understand  the  limiting  behavior  of 
calculations  on  lP3's  and  APB's  in  the  companion  paper.111 
4.  The  Ordering  Limits  at  T^O 

In  the  ground  state  under  conditions  when  w<0  and  A^B  clusters 
exist  1(0,1)  in  Table  1]  there  must  be  a phase  transition  from  the  fee 
(disordered)  to  the  A^B  stiucture  as  x increases  from  0 to  0.25.  The 
energy  provides  uo  clue  since  it  is  linear  over  the  entire  range,  (fig. 

2).  The  classical  common  tangent  construction  based  on  the  energy  does 
not  provide  definite  compositions,  for  it  is  tangent  along  the  entire 
composition  range.  For  the  same  reason  chemical  potentials  are  also 
constant  over  this  composition  range  (fig.  2).  The  terminal  composition 
r= 0 is  not  the  phase  limit  for  obviously  a dilute  mixture  of  B in  A need 
not  be  a two-phase  system.  The  answer  is  not  to  be  found  in  the  energy  but 
in  the  entropy  or  the  ground  state  degeneracy  as  it  affects  the  entropy. 

In  this  section  then  we  first  undertake  to  calculate  the  free  energy  at 
finite  temperature  and  examine  the  low-temperature  limit  under  the 
condition  where  the  (0,1)  limit  of  Table  1 applies.  Later  in  this  section, 
we  derive  the  combinatorial  equations  using  the  tetrahedron  approx- 
imation for  ground  state  degeneracy  and  examine  the  entropy  for  two- 
phase  behavior.  The  two  procedures  give  equivalent  results. 

The  * wo  species  A and  B are  designated  by  i=l  and  2,  respectively. 

In  the  tetrahedron  approximation  of  CVM,  the  basic  variable  z,.^g  *s 
the  probability  of  finding  atomic  species  i,  j,  k and  £ on  the  four 
vertices  of  a tetrahedron.  We  further  require  that  the  fourth  subscript 


LATTICE  LINE  ACROSS  THE  BOUNOARY 


FIG.  10.  Density  profile  versus  the  boundary  calculated  by  the 


sum  method. 


ucu  auuui  me  somewnat 
bothersome  shape  of  the  3x3  curve  for  low  tempera- 
tures. 

One  other  result  of  the  3x4  cluster  calculation  worth 


J.  Chem.  Phys.,  Vol.  68,  No.  1,  1 January  1978 


£ is  for  the  atom  located  on  the  sublattico  which  is  ; iv lwrcnti.il  1. 

occupied  by  B atoms  in  the  ordered  A^B  structure.  For  the  sake  of 

brevity  we  call  this  sublattice  the  B sublattice  and  the  rest  the  A 

sublattices;  the  three  A sublattices  are  equivalent. 

Besides  z^^'s  ve  use  the  pair  variables  y„  and  together 

with  tfce  point  variables  x.  and  u„.  For  y. . both  atoms  are  on  the 

l £ ij 

A sublattices  while  in  v.„  the  first  and  second  subscripts  indicate  the 

i£  r 

species  on  the  A and  B sublattices,  respectively.  For  the  point 
variables,  is  the  probability  of  finding  an  i^*1  species  on  an 
A subls  -:e  point,  and  u^  is  the  corresponding  quantity  on  the  B 
sublattice.  These  variables  are  connected  by  the  geometrical 


relations: 


y . z z 

yij  k,£  zijk£ 


. y 

V — 7 

vi£  j,k  ijk£ 


_ y 

Xi  j ,k,£  Zijk£ 


U£  i,j ,k  Zijk£ 


The  normali'ation  of  z's  is 


i.  £ 2 

1 i, j ,k,£  zijk£ 


(14) 


ci^Oi  iq K„)  across 


nuic  umi  in  «ao  i»*wv»  **«v  v.»r.  - *•*  » • — 

that  the  boundary  free  energy  is  closely  tied  to  the  multiplane  distribution  function  Pn 
the  boundary. 


Now  we  go  to  the  puzzle  mentioned  at  the  end  o£  Sec.  II.  The  key  to  solve  the  puzzle  Is  the  fact  that  the  probabil- 
ity expression  for  m * 1 planes  , v, In  (A3)  requires  that  the  h function  lies  always  on  the  left  of 
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When  the  grand  potential  C 


C 


E - TS 


(15) 


is  writter.  in  terms  of  z^^'s  and  then  is  minimized  with  respect  to 
z's,  we  obtain  the  following  set  of  equations: 


where 


(16) 


Wijki  H - £ijki  + + + + ^)/8 

Xijki  yij  ^ik  ^jk  vi£  vj£  vkJK 
Xijkje  H xi  xj  xk  uf 


(17) 

(18) 

(19) 


Although  the  order  of  subscripts  are  meaningful  in  z^^,  anc* 

Xijk£’  or<^er  i®  ircsaterial  in  v.  .uf  since  the  energy  parameter 


ijk£ 


Eijk£  ^oes  n°t  depend  on  the  order  of  the  subscripts.  The  quantity 
exp  (\/2kT)  in  (16)  is  the  normalization  factor  and  is  determined 
from  the  normalization  condition  (14). 

When  the  energy  parameters  and  the  chemical  potentials 

together  with  the  temperature  T are  given,  the  equilibrium  state 
of  the  system  is  solved  by  finding  z^^'s  which  satisfy  eqs.  (10) 
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through  (15)  simultaneously.  When  these  equations  are  solved,  tin- 
parameter  X is  identified  as  the  grand  potential  G per  lattice  point: 


X = G/K  (20) 

where  X is  the  total  number  of  lattice  points  in  a system. 

lierause  the  system  we  are  calculating  has  a fixed  number  of  lattice 
points  X,  individual  chemical  potentials  are  not  defined;  only  their 
difference  is.  We  can  arbitrarily  choose 


= - “2  ~ - **  ■ 


At  T=0,  the  entropy  part  does  not  contribute  and  we  can  write 

_ 3E/N 

^2  ' 


The  quantity  2p  is  the  diffusion  potential  of  I.arche  and  Cahn. 

It  is  tabulated  in  Table  1 for  the  various  ground  states.  From  Table 
1 and  figure  2 we  see  that  p is  a step  function  of  composition.  In 
the  (0-1)  range  applying  equation  to  the  energy  listed  in  Table  1 we 
obtain 


p = 2EX  = - 6(l+a)  J w | 


Note  that  v<G.  This  must  also  be  the  value  of  p for  the  two-phase 

i 

equilibrium  at  T-0.  At  low  temperature,  we  expand  p away  from  the 
value  in  ,.'d)  and  write 


|J  = - 6(1  Hi)  |v|  ^ ukT  + 

in  which  a is  the  expansion  coefficient  yet  undetermined.  In  the 


(24) 


quantities  w^^  in  (17)  we  use  the  definition  of  the  energy 
parameter  E^^  in  eq.  (4.2)  of  Heforence  7.  The  explicit  form  of 
Wijk£  S are  in  e<l*  (31)  below  in  the  following  section.  For 

the  derivation  of  this  section,  we  use  the  first  two:  and 

W1112’  Since  the  numerator  of  kT  in  eq.  (16)  should  vanish  at  T-0  we 

expand  X also  as 


A = - 6(l+cf)  | w | + bkT  + . . 


Wien  We  use  the  expansions  (24)  and  (25),  we  see  that  w^^  and 
W1112  are  different  by  a/4.  Then  the  general  formulas  in  (16-19) 
reduce  to  the  following: 


21111  “ exp 


(b  a\  (ylJ 

\2  ' 27/  3 
' / ( x,  u 


(hi  h.)3/2 


W “O' 


(h  a\(Vll  v12) 

L2  2 " 5/  / 3 \5/8 

\ /(x1  »2) 


*1121  ~ *1211  z2111  = cxp 


Ml 


( 22  \ 1/2 

vii  yia  vn  n) 

(X1  X2  V 


The  rest  at  z^^^'s  are  negligible  and  are  not  needed.  These 

are  for  tLe  A_B  phase.  The  z...  's  for  the  disordered  phase 
3 r ljki  * 

can  be  cVtained  from  (26)  by  imposing  the  disordering  conditions: 


y . . = v d x.  = u. . 

ij  /.i  i i 


The  formulation  is  now  complete.  We  use  the  computer  t < ■ sol\.  i In- 
set of  equations  in  (26)  for  different  values  of  a.  For  the  compulation 
we  used  the  values  of  o-O.Ol  and  p=-0.0S  which  can  make  the  Cu^Au  phase 
diagram  best  fit  with  experiments  as  will  be  discussed  in  the  accompany- 
ing poper.^  The  quantity  b is  determined  as  a function  of  the  quantity 
a from  the  normalization  of  z's: 


1 = z + z + 3z 

1111  1112  JZ1 121 


Note  that  the  rest  of  Zjjyj's  are  negligibly  small.  The  solid  curve  in 
figure  4 shows  the  result.  We  repeat  the  solution  for  the  disordered 
phase,  and  obtain  the  dashed  curve  in  figure  4. 

Since  A is  the  grand  potential  as  was  mentioned  in  (20),  the  point 
at  which  the  two  curves  cross  in  figure  4 represent  the  coexistance  of 
the  two  phases.  The  value  of  a at  the  intersection  is  a=0.8109  and  is 
the  right  value  of  the  gradient  of  the  kT  vs.  p curve  at  T=0  in  figure 


The  disorder-Cu^Au  phase  boundary  points  at.  T=0  in  figure  6 were 
calculated  in  this  way.  Since  the  exponential  factors  in  (26)  do  noL 
depend  on  T,  the  phase  boundary  curves  in  this  figure  are  vertical  near 
T=0.  It  ij  believed  that  this  vertical  property  is  a general  result. 

Another  general  property  we  can  deduce  from  (26)  is  the  fact  that 
the  equations  are  independent  of  the  parameters  Of  and  P within  the 
limits  imposed  by  (0-1)  in  Table  I.  Because  of  this,  the  positions  of 
t lie  phasf-  boundaries  at  T > 0 are  independent  of  a and  p.  This  is 
rcadil  _ ierstood  since  this  problem  involves  only  different  ways  of 
arranging  r..v.e  same  ground  state  clusters. 


Similar  calculations  were  done  for  the  Cu„Au-CuAn  boundary  jt 

a 

T -*  0 and  the  results  are  plotted  in  figure  6. 

As  was  mentioned  in  the  introductory  paragraph  of  this  section, 
to  start  from  eq.  (16)  is  not  the  only  way  of  deriving  the  T=0  phase 
boundary.  An  alternative  method  is  to  work  with  the  entropy.  For 
the  (0-1)  cluster  pair,  the  ground  state  energy  is  linear  in  the  com- 
position x as  is  shown  in  the  first  row  of  Table  1,  and  thus  does  not 
contribute  to  the  phase  separation.  Therefore,  the  phase  boundaries  at 
T=0  between  the  disordered  fee  phase  and  the  A^B  phase  are  to  be  calculated 
by  the  common  tangent  construction  based  on  the  two  entropy  curves  as 
shown  schematically  in  figure  7. 

We  now  show  that  to  draw  and  as  functions  of  the  composition 
x and  then  determine  the  phase  boundaries  and  x^  from  the  common 
tangent  is  equivalent  to  the  method  presented  in  this  section. 

For  this  purpose,  we  introduce  a parameter  a and  define  functions 
<t>D(a)  and  <t>Q(a)  as 


<t>D(a)  = Max  ^Sp(x)  - axj 
4>0(a)  = Max  |sq(x)  - axj 


(28) 


When  ws  pint  the  two  functions  ^(a)  and  ^(a)  against  the  parameter  a, 
they  cross  as  in  figure  8.  The  point  P at  which  the  two  curves  cross 
corresponds  to  the  common  tangent  situation  in  figure  7,  The  proof  is 
the  following: 


At  I.  we  have 


W = SD(xD)  ' a0XD  = S0(X0)  " Vo  = W (29;,J 


and  since  the  <t>'s  are  maxima 


dSD(xD) 


= a0  = 


d W 


(29b) 


Thus,  we  obtain 


S0(x0)  ~ W 

X0  " XD 


d SD(XD)  d S0(x0) 


This  shows  that  the  two  curves  Sjj(x^)  and  Sq(xq)  posses  a common  tangent. 
Therefore,  x^  and  are  the  boundaries  of  the  two  phases. 

When  we  formulate  the  entropy  using  the  tetrahedron  approximation 
of  the  CVM,  and  proceed  formulating  the  functions  4*^  and  4>q,  we  arrive 
exactly  at  the  eqs.  in  (26). 

One  conclusion  which  is  clearly  derived  from  the  illustration  in 
figures  7 and  8 is  that  the  phase  boundaries  are  determined  by  the 
entropy  formula  only,  and  thus  depend  on  the  approximation  used  in  the 


The  singularity  at  T=0  is  illustrated  in  figures  5 and  6.  Each  of 
the  six  lines  numbered  1 through  6 in  the  two  figures  are  two  represen- 
tations cf  an  approach  to  T=0.  For  any  x in  the  range  in  the  range 
0<x<l/i,  p approaches  the  same  limit,  each  x in  figure  6 corresponding 
to  a fixed  limiting  slope  in  figure  5. 


5.  Altjrr.a~.ive  Derivation  of  Figure  1 

Cocbir.3i.ions  of  figures  5 and  6 and  of  figures  7 and  8 indicate 
the  dual  properties  between  the  p and  composition  space  analysis.  In 
the  present  section  we  briefly  show  an  alternative  derivation  of  figure 
1 based  on  the  p space  analysis. 

We  go  back  to  eq.  (16)  for  In  this  expression  the  Y and 

X parts  ccne  from  the  entropy  expression  and  hence  we  can  disregard 
in  our  derivation  of  figure  1.  The  energy  and  the  chemical  potential 
iuforaiation  is  contained  in  the  quantity  w^^.  Because  of  the 
permutation  symmetry  among  the  subscripts,  there  are  five 
which  are  written  explicitly  here: 


wnn  = " 2 M 

w1112  = ' 2 (1+“)w  " 4 


W1222  = “ 2 (1+P)w  + 5 M 
W2222  = 2 M 

Each  of  the  five  w^jjc£,s  in  (31)  is  a line  in  the  w^^£  vs-  M 
space.  Of  these,  w^  WH22’  an<*  w2222  are  indePen<ient  °f  a and 
p and  are  drawn  by  thick  lines  in  figure  9.  For  the  purpose  of 
illustrate  :r. , figure  9 is  drawn  for  the  case  w<0. 

Ir.  i-j-re  9 for  the  region  |j<-4 1 w | , w^^  is  larger  than  w^.^. 
This  near:  that  in  this  region  of  p,  the  cluster  1111  is  more  stable 
than  the  .luster  1122  in  T=0  because  w„^£  is  divided  by  T.  By  draw- 


: 

II 

iug  five  °f  (31)  on  this  diagram  ajid  comparing  which  is  the 

largest,  we  see  which  cluster  is  present  in  the  system  for  a given  value 
of  p.  When  two  w^^  lines  cross,  it  means  the  two  clusters  can  coexist. 

The  Wjjj2  (31)  depends  on  a and  is  drawn  by  a chain  line  of 
negative  slope  in  figure  9.  The  v±222’  wkich  depends  on  p,  is  drawn  by 
a chain  line' of  positive  slope. 

When  P<-l/3,  w^222  *s  alwaYs  suppressed,  as  indicated  by  the  w^0<>2 
line  marked  by  8^.  This  corresponds  to  the  fact  that  AB^  does  not 
appear  in  the  region  P<-l/3  in  figure  1(a). 

At  the  point  P in  figure  9,  three  lines  meet:  wiu2^a^’  w112">  anc* 
w1222^2^‘  ^or  ^^2  wheD  “ is  in  the  ran8e  -l/3<a<a2,  all  five  clusters 
can  be  stable  at  some  value  of  \i.  The  combination  p=p2  3 11 3 tf=cr 2 is  a 
boundary  such  that  either  a2<d  or  P2<P  makes  the  w^j2->  phase  suppressed. 
At  P,  the  equations 


W1112  “ W1122  W1222  Co“} 

hold  and  use  of  (31)  leads  to  the  relation 

«+P  = f (33) 

which  is  the  line  marked  (P)  in  figure  1(a). 

At  Q in  figure  9,  three  lines  meet:  wuj2^0,P»  w]222^2^  and 
w?.22")'  ^oceedins  in  the  similar  way,  as  in  the  preceding  para- 
graph, we  arrive  at  \ 

t 

a-3p=2  (34) 

which  ii  the  line  marked  (Q)  in  figure  1(a).  The  line  marked  (U) 
in  figure  1(a)  corresponds  to  the  point  R in  figure  9. 


: 


6. 


Summary  and  Concluding  Remarks 

Multiatom  forces  (characterized  by  two  parameters  a and  p), 
rather  than  the  composition-dependent  energy  parameters,  were  used 
recently  successfully  in  deriving  the  asymmetry  of  Cu^Au  phase  diagram.®’^ 
This  formulation  forms  the  basis  of  the  accompanying  paper  ou  phase 
boundaries.1®  As  a study  to  back  up  the  use  of  the  parameters  or  and  p, 
their  effect  on  the  stability  of  phases  at  T=0  is  worked  out  in  the 

first  part  of  the  present  paper.  The  results  are  shown  in  figure  1. 
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In  figure  1,  as  well  as  in  previous  studies  ’ of  possible  phases 
at  T=0,  only  the  energy  expression  of  the  system  comes  into  play.  In 
the  second  part  of  the  paper,  we  derive  the  position  of  phase  boundaries 
at  T-*0.  This  is  done  using  the  entropy  expression,  and  hence  the  position 
of  the  boundary  depends  on  the  approximation  used  in  the  entropy  ex- 
pression. Figures  4,  5 and  6 show  the  results.  A general  conclusion  is 
that  in  the  temperature  T vs.  composition  x plot,  a phase  boundary  near 
T=0  is  always  parallel  to  the  T axis  and  comes  straight  down  to  the  x 
axis. 

The  dual  nature  between  the  composition  and  the  chemical  potential 
p is  one  theme  repeated  in  the  paper.  The  fact  that  p takes  the  same 
constant  value  for  a range  of  x in  the  single  phase  regions  is  illustrated 
in  figures  5 and  6. 
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Table  1 


Clusters 
present  (n,m) 

Range 
in  x 

Energy 

(see  Table  2) 

Chem . 

Pot.  Cm) 

Conditions 
(see  Table 

_2) 

0,1 

0-1/4 

4xE^ 

2E1 

| 012|  >0 

| 013]  >0 

| 014!  >0 

0,2 

0-1/2 

2xE2 

El 

| 012| <0 

| 023|  >0 

| 024] >0 

°.3 

0-3/4 

§**3 

h 

3 3 

| 013|  <0 

| 023|  <0 

| 034]  >0 

0,4 

0-1 

0 

0 

| 014]  <0 

| 024]  <0 

| 034] <0 

1,2 

1/4-1/2 

(2-4x)E1+(4x-1)E2 

2(E2-E1) 

| 012] >0 

| 123| >0 

| 124] >0 

1,3 

1/4-3/4 

(|  -2x)E1+(2x-  |)E3 

| 013] >0 

| 123| <0 

| 134]  >0 

1,4 

1/4-1 

3(1-x)Ej. 

- -E 

3 1 

| 014| >0 

| 124| <0 

| 134] <0 

2,3 

1/2-3/4 

(3-4x)E2+(4x-2)E3 

2(E3-E2) 

| 023] >0 

| 123 | >0 

| 234]  >0 

2,4 

1/2-1 

2(1-x)E2 

-E2 

] 024] >0 

| 124] >0 

| 234] <0 

3,4 

3/4-1 

4(1-x)E3 

- 2E3 

| 034] >0 

| 134 | >0 

| 234] >0 

(a)  is 


* 


C 1 £ . _ 


Fic.  5. 


Fit.  t. 


Figure  Cap t ion s 

Classification  of  stable  phases  in  tin-  n .uni  s 
for  w<0  and  (b)  is  for  w>0. 

An  example  of  the  ground  state  enei gy-cemposit  Lou  inagran. 

This  is  for  the  case  where  A^B  and  AB  are  the  only  intermediate 
phases.  In  the  range  0<x<l/2  solid  solutions  are  passible, 
the  question  of  two  phase  regions  shown  dashed  is  examined  in 
Section  4.  Where  there  is  a missing  intermediate  phase  in  the 
range  1/2<x<1  no  solid  solutions  are  possible  in  the  ground 
state. 

The  chemical  potential  corresponding  to  the  situation  in  fig. 

2.  Note  that  this  is  a step  function.  The  chemical  potentials 
at  two-phase  equilibrium  are  fixed  in  this  construction, 
although  the  coexisting  compositions  are  not. 

Plots  of  b vs.  a from  computer  calculations.  Note  that  the 
disorder  and  the  A^B  phases  cross  at  a = 0.8109.  The  tetrahedral 
cultiatomic  interaction  parameters  chosen  are  a-O.Cl  and  £=- 

o.cs. 

The  chemical  potential  p and  temperature  diagram  calculated 
for  w<0,  a=0.01  and  p=-0.08.  Note  that  p=-(6+u)  v at  T=0. 

The  lines  indicated  by  1 through  6 all  radiate  freer  this 
print. 

Tr.e  composition  and  temperature  phase  diagram  corresponding 
to  figure  5.  The  short  numbered  linos  correspond  to  the 
lines  with  the  same  number  ,in  figure  5. 


1 


'4 


Fig. 


Fig.  3. 


Fig.  9 


Schematic  diagram  of  entropy  curves,  for  the  disordered 
phase  and  Sy  for  the  ordered  A^B  phase.  The  common  tang"nt 
determines  the  phase  x..  and  .\,. 

Schematic  plot  of  the  functions  <t>^(a)  and  <hy(a)  in  eq.  (28). 
The  point  P at  which  the  two  curves  cross  correspond  to  the 
coexistence  of  two  phases. 

Working  diagram  which  is  used  to  lead  to  fig.  1(a)  from  the 
alternative  treatment  of  Section  5.  Points  P,  Q and  R corre- 
sponds to  the  lines  (P) , (Q)  and  (R)  in  fig.  1. 
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Abstract 

Equilibrium  antiphase  (APB)  and  interphase  (IPB)  boundaries  in 
the  Cu-Au  system  are  examined  theoretically  using  the  cluster  variation 
method  with  multi-atom  interactions  whose  magnitudes  were  previously 
obtained  from  a fit  of  the  phase  diagram.  The  IPB  energy  between 
equilibrium  disordered  fee  and  ordered  Cu^Au  phases  is  strongly 
tenperature  (and  hence  composition)  dependent,  being  much  higher  for  the 
copper-rich  side  of  the  congruent  point  T , and  having  a maximum  on 
that  side  at  about  0.95  Tc»  The  IPB  energies  are  only  slightly  anisotropic. 
The  APB  energies  at  constant  chemical  potential  decrease  monotonically 
with  increasing  temperature;  at  constant  nonstoichicmetric  composition 
they  increase  at  lew  temperature  to  a maximum  well  below  the  disordering 
temperature.  Near  the  congruent  point,  the  APB  undergoes  one  or  more 
second  order  surface  phase  transition  in  which  an  interfacial  layer 
resembling  inhomogeneous  Cu-Au  (L1q)  phase  develops  within  the  boundary. 

The  APB  with  (hkO)  orientation  (in  our  notation)  should  be  perfectly  wet 
by  the  disordered  phase  at  the  disordering  temperature  for  that  particular 
composition. 


[ 
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I.  Introduction 

In  ordering  systems  there  occur  two  kinds  of  particularly  simple 
coherent  interfaces,  antiphase  dcnain  boundaries  (APB)  and  coherent 
interphase  boundaries  (IPB).  If  we  ignored  the  difference  in  the 
identity  of  the  atomic  species,  the  two  domains  or  phases  meeting  at 
coherent  interface  would  be  part  of  the  same  single  crystal  as  is 
illustrated  in  Fig.  1.  An  APB  separates  two  dcnains  of  the  same  ordered 
phase.  An  IPB  separates  two  different  phases.  IPB's  can  occur  between 
disordered  phases  differing  in  composition,  between  an  ordered  and  a 
disordered  phase,  or  between  two  ordered  phases. 

APB's  can  exist  in  all  ordered  phases.  They  result  from  the  symmetry 
breaking  during  ordering  processes  which  can  start  in  different  ways  in 
various  locations  in  a disordered  lattice.  The  APB's  form  wherever  two 
such  regions  contact . APB's  result  also  from  the  presence  (or  motion) 
of  dislocations  whose  Burger's  vector  are  not  translation  vectors  of 
the  super lattice.  The  mechanical  properties  of  superlattices  are  strongly 
effected  by  the  fact  that  motion  of  such  dislocations  changes  the  area 
of  APB's.  Often  deformation  can  only  occur  by  groups  of  dislocations 
whose  Burger's  vectors  sum  to  a super lattice  translation  vector,  and 
whose  motion  as  a group  resotres  long-range  order  in  the  structure  [1,2] . 

The  IPB's  are  important  in  alloys  which  order  by  first-order  transition 
The  free  energy  of  an  IPB  is  a factor  in  determining  the  rate  of 
nucleation  [3]  of  the  ordered  phase  on  cooling  or  the  disordered  phase 
on  heating.  For  multiphase  alloys  it  also  affects  the  shape  and  dispersion 
of  particles  and  is  the  main  driving  force  in  the  long-time  coarsening 
of  such  a dispersion  [4]. 


Phase  transitions  within  such  interfaces  have  been  predicted  [5], 
but  the  order  of  the  transition  is  very  high,  so  that  all  properties 


that  are  related  to  low  derivatives  of  the  surface  free  energy  will 
be  continuous  through  this  transition.  There  would  be  important 
implication  of  a change  in  character  of  these  interfaces  on  properties 
of  such  alloys,  if  the  order  of  the  surface  transition  were  lower. 

Whether  the  ordering  transition  itself  is  first  order  or  higher 
order  is  an  important  factor  in  discussing  such  interfaces.  At  the  critical 
temperature  of  a higher-order  transition  the  disordered  and  the  various 
domains  of  the  ordered  phase  all  become  identical  to  each  other. 

Consequently  APB's  disappear  as  the  critical  temperature  is  approached 
and  their  energy  vanishes  [6-10],  Because  there  is  no  coexistence  of 
ordered  and  disordered  phases  (except  trivially  at  the  critical  temperature) 
there  are  no  equilibriun  IPB's  between  phases  related  by  higher-order 
transitions . 

For  first-order  transitions  two-phase  equilibrium  occurs  over  a 
range  of  temperatures  and  the  two  phases  are  never  identical.  The  ordered 
phase  retains  a finite  amount  of  order  up  to  the  transition  temperature. 
Consequently  both  APB's  and  IPB's  exist  right  up  to  the  transition  temperature 
and  there  is  no  reason  to  assume  that  their  energy  would  vanish  there. 

Several  ordering  transitions  are  first  carder.  In  this  paper  we 
shall  examine  IPB's  and  APB's  in  an  fee  (Cu-Au  type)  system  which  undergoes 
an  ordering  transition  to  fcrm  a phase  with  the  Cu^Au  structure  (Llj). 

In  addition  to  the  Cu-Au  system  these  interfaces  also  occur  in  a number  of 
important  nickel-base  systems,  notably  in  the  Ni-Al  system  which  is  the 
basis  of  superalloys.  These  alloys  are  treated  to  form  a fine  scale 
dispersion  of  coherent  Ni^Al  in  a disordered  nickel  solid  solution. 


They  are  useful  because  their  strength  increases  with  increasing  temperatures 
to  approximately  800°C,  and  then  declines  at  higher  temperatures  [11,12]. 
.Increases  of  strength  with  increasing  temperature  seemr  to  be  quite 
commonly  observed  in  ordered  alloys  [13-18]. 

APB's  have  been  studied  theoretically  for  the  CsCl  (or  CuZn) 
structure  derived  from  the  ordering  of  bcc  [8-10].  In  the  model  chosen, 
(near-neighbor  pair-wise  interaction  energies)  this  is  a second-order 
•transition.  The  corresponding  model  for  fee  gave  completely  unrealistic 
phase-diagrams  [19]  until  it  was  recently  shown  that  these  resulted  from 
the  Bragg- Wi  1 1 iams  and  pair  approximation  [18-21]  rather  than  the  model 
itself.  A cluster-variation  method  (CVM)  in  the  tetrahedron  approximation 
can  give  a phase  diagram  in  which  there  are  three  ordered  phases,  A^B,  AB 
and  ABj,  which  disorder  by  first-order  transitions  [20-21].  A close 
natch  to  the  phase  diagram  found  for  the  Cu-Au  system  is  obtained  if 
four-body  interactions  are  introduced  [22].  The  phase  diagram  calculation 
nade  it  feasible  to  undertake  a theoretical  study  of  APB's  and  IPB's  in 
such  a model  system. 

A previous  study  of  IPB's  in  ordering  systems  [23]  suffers  from 
several  assumptions  and  possibly  some  computational  errors  that  have 
limited  its  applicability.  In  this  study  IPB's  were  created  theoretically 
between  phases  of  arbitrary  composition  and  order,  not  necessarily  ones 
that  would  ever  be  in  equilibrium  with  each  ot hereby  cutting  and  joining- 
and  oompxiting  the  energy  from  the  changes  in  the  bond  count.  No 
rearrangement  of  atoms  near  the  boundary  was  permitted.  Such  "weld" 
interfaces  can  have  negative  energies  as  well  as  positive  and  the  energ;  ’ 
can  almost  always  be  reduced  further;  Under  some  conditions  rearrangement  might 
tend  even  to  complete  homogenization  by  rearrangement.  We  note  that 
for  many  of  the  ground  state  cases  we  obtain  rigorous  results  that  differ 
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frcm  those  obtained  previously  [23]  even  When  we  assure  the  same  interaction 
energies.  The  study  also  examined  APB  energies  without  regard  to  whether  or 
•not  the  initial  phase  oould  be  in  equilibriun.  For  the  Llj  phase  negative 
APB  energies  were  reported  [23]  under  condition  where  D022  should  have 
been  the  stable  ground  state  structure  [24].  Here  the  ”APB"  was  the  first 
- steo  towards  equilibriun. 

In  this  paper  we  will  use  the  model  with  four-body  interaction  energies 
chosen  to  natch  the  Cu-Ab  phase  diagram,  and  calculate  properties  of 
APB’s  and  IPB’s  which  have  been  allowed  to  equilibrate  with  respect  to 
an  atomic  rearrangement  in  a crystal  constrained  to  have  two  different 
domains  or  phases.  At  high  temperatures  and/or  slow  defamation  rates 
there  is  time  for  diffusional  rearrangement.  Such  equilibrated  boundaries 
obey  all  the  laws  of  surface  thermodynamics  including  the  Gibbs  adsorption 
equation  which  can  be  used  for  a sensitive  check  on  the  calculation 
procedures. 

If  an  APB  is  created  so  rapidly  by  shearing  that  there  is  not  time 
far  diffusional  rearrargement  of  atoms,  we  obtain  a very  different 
boundary.  Calculations  of  the  properties  of  such  boundaries  have  been 
made  for  various  ncdels  [7,25]  and  will  be  examined  in  a subsequent 

‘ 

paper  for  the  model  used  here  [26]. 

I 


IX. 


1 


Crystallography  and  Definition  of  Surface  Orientation  Indices 


Tha  disordered  phase  is  foe  (Iln3m)  with  four  equivalent  sites  with 
m3m  point  sysnetry  at  000,  1/2  1/2  0,  1/2  0 1/2  and  0 1/2  1/2.  The 
ordered  phase  has  the  Cu^Au  crystal  structure  (ftn3m)  with  one  kind  of 
occupation  for  the  site  arbitrarily  chosen  at  000  with  m3m  syrtnetry  and 


another  occupation  at  the  three  equivalent  sites  with  4/nmn  symmetry  at 
1/2  1/2  0,  1/2  0 1/2  and  0 1/2  1/2  [27],  Four  domains  are  possible 
because  any  of  the  four  equivalent  sites  in  fee  could  have  become  the 
origin  far  the  Cu^Au  structure. 

It  has  been  customary  in  this  field  to  describe  the  fee  structure  by 
four  inter-penetrating  simple  cubic  "sublattices"  each  centered  on  one 


of  the  four  equivalent  sites  (Fig.  2).  The  Cu^Au  type  ordering  occurs 
when  the  occupation  of  one  of  the  sublattices  becomes  different  from  the 
occupation  on  the  other  three. 

We  define  our  crystallographic  axes  along  the  cube  axes  common  to  both 


structures.  Far  the  IPB  we  choose  the  origin  of  the  coordinate  system 
such  that  it  coincides  with  the  atom  having  m3m  point  symmetry  in  the  ordered 


phase.  For  the  APB  we  place  the  origin  at  this  point  in  one  of  the  domains 

and  arbitrarily  choose  the  z-axis  so  that  the  m3m  position  occurs  at 

1/2  1/2  0 in  the  other  domain  (See  Fig.  3).  The  orientation  of  the  boundary 

A 

is  specified  by  its  normal  n,  but  all  such  properties  must  be  invariant 
to  symmetry  operations  of  a point  group  which  is  cannon  to  both  domains 
or  phases.  For  the  IPB  thip  cannon  point  group  is  m3m.  For  the  APB  it  is 
4/nmn  because  the  translation  which  carries  the  occupation  from  the  site 
at  000  to  the  1/2  1/2  0 site  destroys  the  three-fold  axes  and  two  of  the 
fbur-fold  axes  (Fig.  3),  Accordingly  while  (100),  (010)  and  (001)  IPB's 


I 


; 


are  completely  equivalent,  (0Q1)  APB's  differ  from  (010)  and  ilOO). 

The  former  is  called  a conservative  APB  because  it  could  hove  been 

created  from  a single  domain  by  conservative  slip;  while  the  (100)  and 

(010)  ore  nonconservative,  and  with  perfect  order  either  two  enriched  or 

depleted  layers  are  adjacent  at  the  APB.  They  have  also  been  called 
e 

bpndaries  of  the  first  and  second  kind.  We  may  represent  all  scalar 
surface  properties  by  plots  on  the  unit  sphere,  but  the  symmetry  makes 
the  triangle  with  corners  at  (001),  (101)  and  (111)  sufficient  to  describe 
all  orientations  for  the  IPB,  while  the  larger  triangle  with  earners  at 
(001),  (100)  and  (110)  is  necessary  for  the  APB's. 

The  symmetry  properties  of  bicrystals  with  planar  boundaries  have  recently 
been  developed  [28].  The  formulation  is  general  enough  to  include  APB's 
(no  orientation  change,  just  a translation  of  one  crystal  relative  to  the 
other)  and  IPB's.  Application  of  this  new  group- theoretical  method  to  the 
present  system  leads  to  equivalent  descriptions. 


III.  General  Formulation 

In  this  section  we  develop  the  formulationa  of  the  (100)  APB  as  an 
example.  The  formulation  for  other  orientations  differ  only  in  the  interactions 
and  configurations  of  atomic  planes  parallel  to  the  interface.  The  APB 
can  exist  over  the  entire  range  of  temperature  and  composition  (or  chemical 
potential ) in  which  the  ordered  phase  is  stable.  For  IPB's  there  is 
the  limitation  on  chemical  potential  imposed  by  the  phase  rule  that  different 
phases  must  coexist  at  the  interface;  For  each  temperature  where  the 
phases  coexist,  an  IPB  occurs  at  discrete  chemical  potentials. 

We  shall  call  the  four  simple  cubic  sublattices  I,  II,  III  and  IV  as 
in  Figure  2.  We  shall  call  the  AjB  domains  in  which  the  I sublattice  is 
preferentially  occupied  by  B atoms  the  domain  I ; the  other  domains  are 


correspondingly  defined.  The  [001]  axis  is  then  determined  by  our  convention. 
The  (100)  APB  Is  formed  when  domains  I and  II  are  pieced  on  opposite  sides 
of  the  boundary  and  allowed  to  relax  to  the  equilibrium  configuration. 

By  the  definition  of  dona  ins  I and  II,  the  III  and  IV  sublattices 
are  preferentially  occupied  by  A atoms  in  both  phases  while  the  occupancy 
pattern  changes  between  the  I and  II  sublattices  as  we  go  from  domain  I 
to  II.  Therefore  we  see  that  the  entire  system  can  be  regarded  as  a 
layer  structure  parfpjendicular  to  the  boundary  consisting  of  P layers 
mostly  of  A atcms  and  Q layers  composed  of  mixed  A and  B atoms. 

We  call  lattice  planes  parallel  to  "the  boundary  as  "parallel"  planes, 
for  short;  they  are  numbered  . . . , n,  n+1,  . . . , as  in  Fig.  2. 

In  formulating  the  free  energy  including  the  boundary  layer  using 
the  cluster  variation  method  (CVM),  we  take  a tetrahedron  like  the  I-II-III-IV 
connected  in  Fig.  2 as  the  basic  cluster.  Each  tetrahedron  is  made  of 


two  points  on  one  parallel  plane  (like  I and  III  on  the  n**1  plane)  and 
of  t*o  points  on  an  adjacent  parallel  plane  (like  II  arvl  IV  on  the 
(n+l)th  plane).  Also  we  note  that  one  of  the  two  points  (like  III)  on  one 
parallel  plane  lies  on  a P layer  and  the  other  point  (like  I)  lies  on  a 
Q layer. 

In  a binary  system  an  A atom  and  a B atom  are  called  the  1st  and  the 
2nd  species,  respectively,  for  mathematical  convenience.  We  let  i^  denote 
the  ith  species,  (i=l  or  2)  on  a P layer  point  on  the  plane.  Suppose 
the  HI,  I,  IV  and  II  points  connected  as  a tetrahedron  in  Fig.  2 are 
occupied  by  the  i*\  k**1  and  l^1  species,  respectively.  Then  using 

our  notation  we  can  say  that  this  connected  tetrahedron  has  the  configuration 
*Ph  ^Qn  ^(n+l)  ^(n+l)*  We  n®  introduce  the  basic  variable  zn  (i,j,k,A)  of  the 
theory;  this  variable  designates  the  probability  that  the  configuration 
appears  in  the  connected  tetrahedron  of  Fig.  2. 

For  the  pair  variable  there  are  five  kinds  as  defined  in  Table  1. 

Note  that  the  crder  of  argunents'  has  significance.  Table  2 defines 
the  probability  variables  for  a lattice  point. 

There  are  geometrical  relations  among  the  variables: 

y (i,j)  = E z (i,j,k,l) 

k,l  n 

VPP*/1*^  = ? >fc»f) 

(1) 

Vpq^ijl)  = E zn(i,j,k,t) 

vQPr/^*^  = ? zn^i»j»k,f) 
i,t 


vQQn<1,l)  * zn(i»j*k»t) 

.x^i)  - E *n<i,j,k,l)  (2) 

j»k»t 

Also  we  have  continuity  constraints 

y_(i,j)  = E z(i,j,k,t)  = E z (k,t,i,j>  (3) 

n k,t  n k,t  n A 

and  the  normalization  of  z is 

1 = E z <i,j,k,l>  (4) 

‘ i»j»k»l  n 

Me  have  thus  defined  the  variables . Now  we  go  on  to  write  the  free 
energy  of  the  system.  The  number  of  lattice  points  in  a plane  parallel 
to  the  boundary  is  written  as  N.  Then  the  lumber  of  tetrahedra  connecting 
the  nth  and  the  (n+l)th  planes  is  2N.  The  total  energy  of  the  system  is 
then  written  as 


E = 2N  £ E e(i,j,k,l)  zn(i,j,k,l) 
n i,j,k,t  n 


where  e(i,j,k,t)  is  the  energy  per  tetrahedron  and  is  written  as 

cd.i.i.D  = o.o 

c(l, 1,1,2)  = 1.5  w(l+a) 

c(l,l,2,2)  = 2.0  w (6) 

c(L,2,2,2)  = 1.5  wCl+5) 
c(2,2,2,2)  = 0.0 

In  these  expressions,  which  are  the  same  as  Eq.  (4.2)  in  KiVuchi-de  Fontaine’s 
paper  [22],  the  order  of  arguments  in  e is  immaterial.  The  parameter  w is 
for  an  interaction  of  an  A-B  pair,  and  its  magnitude  is  related  to  the 


F 

t«perature  scale.  The  two  distent  ionless  parameters  a and  0 represent  the 
four-body  interactions  and  make  the  phase  diagram  asynmetric.  The  values 
chosen  for  the  best  fit  [29]  for  the  observed  Cu-Au  diagram  are 
w/k  = -663°K 

A 

o * 0.01  (7) 

A 

0 = -0.08 


which  lie  well  within  the  range  of  values  that  give  the  proper  ground 


2 

-I  y,N.  a -(N/4)  Z I (y,+  Ml*  yv+  y#)z (i,j,k,l)  Cll) 

i*l  1 1 n i,j,k,l  1 J k 1 n 

Since  we  have  only  two  substitutional  species  in  the  system  and  we  have 
no  vacancies,  we  oan  choose  without  loss  of  generality 
-W1  s V2  5 4 (12) 

where  y is  twice  the  diffusion  potential  defined  by  Larch£  and  Cahn  [31]. 
In  calculating  the  equilibrium  state,  we  fix  T and  y,  and  find  the 

A A 

minimum  of  G [30].  The  minimization  of  G leads  to  the  following  set  of 
equations: 

zn(i,j,k,l)  = z<0)  (i,j,k,l)  exp  Lj  An  + an(i,j)  - an+1(k,t)]  (13a) 

where 

z<°>  (i,j,k,t)  = exp  [-0e(i,j,k,t)  + g-  g(yi+  y..+  uk+  y^)]  x 

t>rn(i’3)vPE»i(i>k)vRJn(i''l)vQfti<5>k)vQqn<5>t)Vl0c>t):|1/2  x CVl) 

Vi>Vi<klVi“>1'S/'  <13b> 

In  (13a),  the  quantity  Aft  is  the  Lagrange  multiplier  for  the  normalization 
in  (4)  and  has  the  meaning,  after  the  iteration  has  converged,  that  the 

A 

grand  potential  G of  the  systan  is  written  as 

G = N Z\  . (14) 

n n 

Note  that  N is  the  number  of  lattice  points  within  one  parallel  plane. 

The  quantities  c*n(i,j)  and  an+1(k,t)  in  (13a)  are  the  Lagrange  multipliers 


(c*n(i»j))  i»  called  the  minor  iteration.  In  this  process  tie  can  choose, 
without  loss  of  generality, 

an(l,l)  = 0 far  every  n (15) 

because  of  the  normalization  on  y's: 

£ y (i,j)  = 1 for  every  n (16) 

i.l  " 

The  details  are  similar  to  those  in  Eq.  (4.14)  of  the  b.c.c.  boundary 
paper  [$]. 

When  a^(i , j ) 's  are  thus  solved  the  output  of  one  major  iteration 

step  is  obtained  from  (13a).  This  output  (zn(i,j,k,l)}  is  used  as  the 

input  for  the  next  major  iteration  step.  This  iterative  solution  of  z's 

has  been  called  the  Natural  Iteration  Method  (NlM). 

As  was  discussed  in  length  before  [22,29,30],  the  NIM  has  the 

* 

advantages  (i)  that  the  grand  potential  G always  decreases  at  each  iteration 
step,  (ii)  hence  that  it  always  converges  whatever  initial  values  the 
iteration  may  start,  and  (iii)  that  the  formulation  avoids  subtraction  so 
that  all  the  variables  «L1  ways  stay  positive  however  small  the  value  may  be. 

For  the  boundary  computations,  we  choose  a system  made  of  either  60 
or  80  parallel  planes.  The  number  N of  atoms  in  one  parallel  plane  drops 
out  when  we  assume  that  N is  sufficiently  large  and  neglect  the  end  effect. 
For  the  initial  conditions,  we  choose  the  left  half  of  the  system  to  be 
in  domain  I and  the  right  half  in  domain  II.  With  this  initial  condition, 
the  interaction  converges  after  about  a thousand  major  iterations, 
sometimes  less  sometimes  more,  and  the  number  of  minor  iterations  is 
less  than  ten  except  near  the  beginning  of  iterations. 

When  the  iteration  has  converged,  we  know  the  structure  of  the 
boundary  expressed  by  the  set  of  the  probability  variables  (zn(i,j  ,k,t )} . 


The  density  profile  across  the  boundary,  for  example,  is  expressed  by 

the  point  variables  (x  (i)}  and  (u  ( j )}  of  Table  2.  Examples  are  shown 

n n 

in  later  sections. 

The  excess  free  energy  a is  defined  [5,10,32}  as  the  excess  of  the 

A 

grand  potential  6 in  the  inhomogeneous  system  containing  the  boundary 

A 

over  the  grand  potential  6^  in  the  homogeneous  phase,  and  o is  normalized 


to  unit  area: 

aNo=G-Gh  = NE  UR  - \J  (17) 

n a 

where  \m  is  the  value  of  in  the  homogeneous  system,  and  a is  the  area 

per  fee  lattice  point  in  the  parallel  plane. 

Derivative  Quantities 

In  a recent  paper  [32],  one  of  the  authors  prepared  a general  formula 
to  estimate  derivatives  of  the  IPB  free  energy.  One  of  the  formulas  is 


,3a. 


[S] 

[N^ 

[N2] 

S’ 

V 

V 

S" 

• Nl” 

V 

e.  — 3 


C\t) 


V 

V 


V 
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In  this  expression,  [S]  is  the  entropy  contained  in  a cylinder  which  is 
perpendicular  to  the  boundary  and  whose  cross-sectional  area  is  unity. 

This  cylinder  extends  far  enough  into  each  phase  that  it  encloses 
sane  of  each  homogeneous  phase  beyond  the  influence  of  the  boundary.  The 
quantity  [N^l  is  the  total  lumber  of  the  jLth  species  (i=l  or  2)  contained  in 
this  cylinder.  The  quantity  S'  is  the  entropy  per  unit  volume  of  the  bulk 
homogeneous , and  N^'  is  the  number  of  the  ith  species  in  a unit  volume  of 
the  bulk  phase  I.  The  [S/N^f^]  is  identified  with  the  surface  excess 
entropy  [32]. 


tte  apply  the  expression  (18)  to  the  (100)  IPB,  far  which  the  pruned 
phase  is  the  disordered  phase  and  the  double- primed  phase  is  the  A^B  ordered 
ptase.  For  this  problem,  we  oan  simplify  (if)  using  the  relation 
Mj'  ♦ N2*  * 1 

and  the  corresponding  equation  for  . 

Since  the  formulation  in  the  present  paper  is  based  on  the  lattice 
structure,  we  oan  write  [S]  and  [N^]  as  suns  over  parallel  planes  n.  By 

i 

using  the  relations  in  (19),  we  oan  reformulate  the  general  expression 


(18)  as 


4 <5T>  * - * tSJn 

n 

where 


(20) 


(21) 


[S]n  5 -(N2*  - N2")_1  CN2*  S"  - N2-  S’  ♦ N2jn(S’  - S")]+Sn 

In  this  expression  Sn  and  N2  n are  the  entropy  and  the  number  of  B atoms 
(i.e.,  the  atom  of  the  2nd  species)  per  lattice  point  on  the  nth  parallel 
plane.  Quantities  S',  S",  N2\  and  N2"  are  per  lattice  point  of  the 
homogeneous  phases . 

A modification  of  this  formulation  oan  be  used  for  APB's  [32],  tut 
because  the  two  homogeneous  domains  have  identical  thermodynamic  properties 
there  is  little  difficulty  defining  invariant  surface  excess  quantities  [10]. 
Ground  State  for  IPB's  and  APB's 


As  in  the  previous  work  on  6-trass  APB's  [8],  it  is  instructive  to 
examine  low  tenperature  variation  of  boundary  properties  with  the  orientation 
and  compositional  or  chemical  potential,  lb  this  we  examine  the  ground  state 
and  ground  state  degeneracy  in  this  model  of  a system  containing  ATB's 
and  IPB's  of  various  orientation.  While  no  real  system  could  reach 
equilibrium  at  low  temperatures,  the  calculated  low  temperature  equilibria 
display  a clue  to  the  understanding  of  trends  that  begin  at  higher  temperatures. 


We  use  the  methods  and  notation  of  the  preceding  paper  [33] . For  the 

A /V 

case  considered  here,  w<0,  a and  6 snail  we  have  for  nonstoichionetric  alloys: 

A A A A 

■ -2(1+3o)wZ2  - (2+3a-$)wZ^  - 12(l4a)wZ^  when  0<x<Jt 

E-(2-4»)E1  - (4x-1)E2  - -2(l+3a)wZQ  - (2-3a-38)wZ3  - 6(l-a)wZA  when  k<x'*h 

where  x is  the  nole  fraction  of  B atoms.  For  the  ground  state  the  r.h.s.  of 
these  equations  is  zero,  because  the  atons  can  be  arranged  with  Zj^Zj-Z^O 
for  the  first  case  and  with  Zq-Zj^Z^- 0 for  the  second.  If  an  IPB  or  APB 
requires  the  presence  of  these  higher  energency  clusters  the  Increase  in  energy 
per  unit  area  is  the  limiting  low  temperature  surface  energy  for  such  a boundary. 

It  is  apparent  that  IPB’s  between  fee  and  A3B  and  between  AgB  and  AB 
can  be  created  without  introducing  any  higher  energy  clusters.  It  is 
readily  demonstrated  by  a construction  such  as  Fig.  1 that  (100)  and 
(001)  APB  can  be  created  entirely  from  allowed  clusters  and  thus  with  zero 
energy.  APB's  with  other  orientation  can  always  achieve  zero  energy  by 
facetting  to  cube  orientations.  Therefore  ground  state  APB  energies  at 
T=0  for  nonstoichionetric  crystal  are  zero  for  all  orientations. 

Stoichiometric  A^B  alloys  can  contain  only  A3B  clusters.  The 
formation  of  (001)  APB  requires  no  energy  because  it  requires  no  other 
clusters,  tut  for  all  other  orientations  other  clusters  are  required.  It 
is  convenient  to  formulate  the  stoichiometric  case  in  terms  of  chemical 
potentials . 

The  result  for  the  (100)  APB  is  sketched  in  Fig.  4 as  a three-dimensional 
stereographic  drawing  of  the  (a,p ,x)  relation  together  with  its  three 
projections.  The  lines  in  Fig.  4 are 

°a2 

w 


fa  _ ( -2(l-3a)  - y/|w[  when  -4  < p/|w|  < -2(l-3a) 
wf  ~ \ P/H“  6(1+0)  when  -6(l+a)  < p/|w|  < -4 


(22) 


and  a is  the  length  of  the  cube  edge  as  shown  in  Fig.  2.  These  relations 
are  derived  by  counting  bends  of  the  APB  between  two  perfectly  ordered 
domains;  we  skip  the  details  of  derivations  and  present  only  their 
interpretations.  The  (y,x)  projection  shows  a step  function  in  accordance 
with  Fig.  3 of  the  preceding  paper  [33],  the  (o ,x)  projection  consequently 
shows  o as  a S -function  of  the  composition  at  the  stoichiometry . The 
(cr,y)  projection  reflects  that  the  adsorption  by  this  APB  is  ±1  B atom  per 
sectional  area  a , since  OCf/3y)T  is  proportional  to  the  excess  B atoms 
adsorbed  on  the  APB.  On  the  w<j-4 1 w | side,  the  boundary  is  composed  of 
a layer  of  A4  tetrahedna,  and  on  the  -4|w|<y  side  there  are  A2B2 
tetrahedra.  At  y=-4|w| , a is  a maximum  and  either  A^  or  A2B2  tetrahedral 
clusters  give  the  same  O’.  The  constant  adsorption  on  the  two  legs  of  the 
(o,y)  projection  are  thus  consistent  with  the  Gibbs  adsorption  equation. 

The  two  values  of  y at  the  bottom  of  the  legs  are  resp.  the  values  at  T=0 
for  the  disorder  - A^B  phase  coexistence  and  far  the  AjB-A^  phase 
coexistence.  The  a converges  to  the  point  (T=0,  x=l/4)  non-unifomly  as  is 
seen  in  Fig.  5,  and  the  nature  of  the  non-uniformity  is  that  the  converged 
values  are  different  depending  on  y as  is  illustrated  near  T=0  in  Fig.  6, 
as  well  as  in  Fig.  4. 

The  energy  of  APB's  of  other  orientation  in  a stoichiometric  alloy  is 

a strong  function  of  orientation.  A polar  plot  of  a vs.  orientation  can  be 

described  [8]  as  a rasberry  figure  composed  of  the  outer  envelope  of  four 

spheres  with  diameters  equal  to  o (100)  and  centers  at  (*—  a (100),  0,  0) 

2 

and  (0,  a (100),  0) 


IV.  Results 


The  results  of  the  calculations  yield  the  profiles  and  thermodynamic 
properties  of  IPB's  and  APB's  as  a function  of  orientation  n in  the 
teRyMrature-conpos  it  ion  (or  temperature-chemical  potential)  region  where 
these  interfaces  could  exist.  APB's  could  exist  over  the  cloned  region 
on  the  phase  diagram  where  the  AjB  phase  is  stable  as  shown  in  Fig.  5; 

IPB's  can  exist  only  at  that  portion  of  the  boundary  of  this  region  where 
both  ordered  and  disordered  phases  are  stable  and  can  coexist . The 
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properties  of  APB's  are  determined  by  three  variables  (n,T,x)  or  (n,T,p) 
since  p»^.(T,w)  Along  the  two-phase  coexistence  curve  p and  the  x's  of  the 
two  phases  are  functions  of  T as  shown  in  Fig.  5 of  the  preceding  paper; 

A 

therefore  IPB  properties  are  determined  by  two  varialbes;  n and  either  p, 

T,  or  the  x's  of  either  phase. 

IPB 

Figures  7 and  8 show  a for  (100)  and  (110)  IPB  as  a function  of  T and 
of  p,  respectively.  It  is  to  be  noted  that  there  is  only  slight  anisotropy 
with  (110)  surfaces  having  slightly  higher  o.  The  anisotropy  falls  well 
below  the  range  that  would  result  in  facetting.  The  equilibrium  shape  of 
an  ordered  domain  within  a disordered  iratrix  would  be  a "sphere''  slightly 
flattened  (by  about  10%)  along  the  (100)  directions.  This  is  consistent  with 
transmission  electron  microscope  observations  that  small  coherent  ordered 
Cu3Au  and  NijAl  in  disordered  matrixes  are  close  to  spherical  [34-36].  As 
the  particles  grow  they  tend  increasingly  towards  rounded  cubes.  Several 
factors  in  addition  to  surface  free  energy,  particularly  coheix?ncy  stmin 
are  important  in  determining  particle  shape,  but  surface  free  energy  is 
dominant  for  small  particles  [4],  It  is  also  to  be  noted  that  neither  a 
minimum  nor  a maximum  in  o occurs  at  the  congruent  (i.e.  maximum) 


i 

twperature  of  the  phase  diagram.  Indeed  because  the  denominator  of  the 
r.h.s.  of  equation  (18)  ia  aero  at  this  maximun  temperature,  do/dT  is 
infinite  unless  the  numerator  is  also  zero  at  such  a point.  Unless  there 
are  special  conditions  the  nunerator  will  not  be  zero  and  such  an  infinite 
slope  at  the  congruent  point  is  a general  result. 

The  value  of  o for  an  IPB  is'  important  in  the  nucleation  of  the 
ordered  phase.  One  prediction  is  that  the  Cu^Au  alloys  enriched  in  Au 
will  nucleate  the  ordered  phase  more  easily,  i.e. , at  less  undercooling 
than  alloys  enriched  in  Cu. 

Consistent  with  our  model  which  has  no  interactions  between  neighbors 
other  than  the  first,  o tends  to  0 at  0 °K.  In  a real  system  low  temperature 
o values  must  depend  on  other  interactions. 

Thermodynamic  self-consistency  within  the  model  is  demonstrated  by 
Figure  9 in  which  o far  (100)  IPB  is  given  at  several  values  of  T together 
with  a slope  given  by  minus  the  surface  excess  entropy  [S/N^N?]  obtained  at 
that  temperature  frem  the  same  IPB  profile  used  to  calculate  o in  Fig.  6. 

This  demonstration  that  equation  (18)  holds  is  a very  sensitive  test  of 
self-consistency  110,37]. 

APB 

Because  of  the  non-uniform  convergence  of  y at  the  point  (T,x)  = 

(0,  1/4),  we  give  in  Fig.  6 o vs  T for  (100)  boundaries  both  at  constant  x 
and  at  constant  y.  At  x = 1/4  and  T = 0,  y can  range  within  the  limits 
given  in  equation  (22).  All  the  constant  y curves  converge  to  x = 1/4  at 
0 °K  and  to  finite  values  of  o,  while  all  the  x / 1/4  curves  lead  to 
o : 0 at  0 °K.  It  is  to  be  noted  that  the  curves  for  a constant  nonstoichiomet rio 
x show  a nuximun  in  o. 

For  the  nonstoichiametric  constant  x curves,  the  major  term  in  the 
low  temperature  dependence  of  o ocmes  from  the  adsorption.  Because  the 


reversible  creation  of  APB  removes  excess  species  from  the  ncnstoichi  erne  trie 
ordered  phase  where  they  have  high  partial  molar  entropy  and  concentrates 
than  alor^  the  boundary  with  low  entropy,  the  net  entropy  change  for  the 
system  is  negative  when  APB  area  is  reversibly  increased.  As  in  the  case 
for  bcc  APB's  [10]  the  maxinun  in  o occurs  in  the  temperature  range  of 
desorption.  Consistency  with  the  Gibbs  adsorption  equation  is  demonstrated 
in  Fig.  10,  to  be  oonpared  with  Fig.  6. 

The  anisotropy  of  APB*  is  pronounced.  In  this  model  the  (001)  APB  has 
o = 0 at  all  temperatures.  The  (100)  and  (110)  results  are  oenpared  in 
Fig.  6.  The  ratio  of  o11Q  to  o100»  and  therefore  on0  does  not  facet  into 
°100  and  °010‘  ®ur  oalculatrion  for  (101)  led  to  a value  of  a which  at  all 
temperatures  exceeded  A/5  o100*  T**  APB  could  always  reduce  their 

energy  by  facetting  to  (001)  and  (100).  Although  we  have  not  computed 
ether  orientations  we  would  expect  the  anisotropy  to  be  such  that  all 
(h,k,t)  orientation  would  facet  to  (h,k,0)  and  (001),  while  (h,k,0)  would 
be  stable. 

The  domain  structures  observed  for  Cu^Au  shew  both  facets  and  rounded 

APB's  that  are  mostly  near  cube  planes  [37].  We  expect  that  the  sharp 

edges  connect  a (100)  with  a (001),  While  the  rounded  portions  lie  along 

(h,k,0),  connecting  (100)  with  (010)  segments.  The  (111)  APB  energies  have 

been  determined  experimented ly  at  350°C  in  a stoichiometric  alloy  to  be 

25  mJ/m2.  We  calculate  aa2/|w1  =0.4  for  this  temperature  far  (h,k,0)  APB's, 

which  is  given  23  rrd/m.  Because  of  the  strong  anisotropy  we  calculate 
2 

18  mJ/m  fer  (111)  orientations. 

Perfect  Wetting  and  Fhase  Changes  in  APB's 

Our  computed  values  of  oa  /|w|  at  the  congruent  point  are  0.1469  and 
0.15,41  for  the  (100)  and  (110)  APB's  respectively  and  0.0722  and  0.0767  for 
the  (100)  and  (110)  IPB's  respectively.  Taking  into  consideration  of  the 
special  limitation  in  the  capability  of  the  computer  in  calculating  the 


APB  congruent  point,  we  oonclude  that  the  canputer  results  indicate 
for  (h,k,0)  orientations  at  the  congruent  point 

°APB  = 2oIPB  (23) 

Thus  the  (100)  APB  could  either  consist  of  two  (100)  IPB's  with  a disordered 
layer  in  between  or  it  coincidentally  has  the  same  free  energy.  Figures 
11a  and  lib  show  the  calculated  profiles  of  the  APB  and  the  IPB,  and 
confirm  that  the  disordered  layer  is  part  of  the  APB  and  that  indeed  the 
APB  has  become  two  IPB's.  There  is  no  reason  why  the  thickness  of  the 
equilibrium  disordered  layer  exactly  at  the  congruent  point  TQ  should  not 
be  infinite;  the  finite  thickness  shown  in  Fig.  11a  is  the  result  of  the 
computer  calculation  which  is  difficult  to  attain  perfectly  converged 
solution  in  the  finite  canputer  space  allocated  far  the  APB  at  Tc- 

We  determined  OT  fcr  (100)  and  (110)  IPB's  from  kT/|w|  = 0.50  on  the 
low  Au  side  of  the  oongruent  point  to  the  eutectoid  temperature  on  the 
high  Au  side  (Fig.  5)  and  compared  this  with  the  limiting  values  at  the 
two-phase  boundary  of  the  corresponding  APB's.  Equation  (23)  held  everywhere. 

The  implication  of  this  result  is  that  in  two-phase  alloys  the  disordered 
phase  will  coat  all  (h,k,0)  APB's  near  T . Since  at  equilibrium  all  other 
orientations  will  facet  into  (h,k,0)  and  (001)  APB’s  there  will  be  only 
(001)  APB's  within  equilibrated  particles  of  the  ordered  phase.  By 
analogy  to  the  corresponding  phenomenon  in  fluids,  this  is  called  perfect 
wetting  of  the  (h,k,0)  APB's  by  the  disordered  phase.  Perfect  wetting  has 
been  shewn  theoretically  to  cocur  near  all  tricritical  points  [39]  and  it  hao 
been  found  that  in  the  Fe-Al  system  the  disordered  phase  coats  all  Fe-Al  APB's 
[HO, 41],  The  general  proof  which  applies  near  tricritical  points  is  not 
applicable  to  the  AjB  APB's  and  we  do  not  know  whether  the  firming  of 
perfect  wetting  is  a general  result. 

Figure  11  shows  the  internal  structure  of  (100)  APB  and  IPB  at  the  congruent 
point.  Within  each  there  is  a thick  layer  in  which  each  plane  parallel  to 


the  boundary  is  disordered  (a  pair  of  sublattices  have  equal  occupations) 
but  adjacent  planes  are  nodulated  in  occupation  by  B.  This  layer  has 
the  CuAu  (LLq)  structure  in  Which  the  CuAu  order  parameter  gradual 
changes  sign  with  distance  perpendicular  to  the  boundary.  At  the  center 
of  the  APB  the  order  parameter  is  zero,  all  four  sublattices  having  equal 
occupation.  The  notation  of  Shockley  [19]  lists  structures  of  phases 
by  the  number  of  sublattices  having  identical  occupation;  (4)  is  disordered 
fee,  (3,1)  is  CUjAu,  (2,2),  is  CuAu.  Using  this  notation  we  find  that 
at  the  congruent  point,  the  layers  of  the  (001)  IPB  in  Fig.  lib  change  from 
(3,1)  on  the  far  right  to  (2,1,1)  to  (2,2)  to  on  the  far  left  (4)  vhile 
the  (100)  APB  consists  of  two  such  IPB’s. 

As  we  lower  the  temperature  at  constant  composition  or  chemical 
potential  the  APB’s  diminish.  The  (100)  APB  profile  is  shewn  for  two 
temperatures  in  Fig.  12a  and  b.  From  the  congruent  temperature  dewn  to 
KT/|w|  = 0.90  the  profile  narrows  but  the  sequence  of  the  layers  remains 
the  same.  At  KT/| w|  = 0.87  the  APB  has  lost  the  (2,2)  and  (4)  layers  and 
contains  (2,1,1)  and  (1,1, 1.1)  layers  instead.  Below  this  temperature, 
there  is  further  narrowing,  but  this  structure  is  retained  to  quite  low 
temperatures. 

In  order  to  ascertain  whether  or  not  there  is  a surface  phase  transition 
between  these  two  temperatures  we  give  in  Fig.  13  the  values  of  the  adsorption 
and  the  excess  in  entropy  far  the  (1001  APB  as  a function  of  temperature . 

There  is  a break  in  the  slopes  of  these  two  curves  at  KT/|w|  = 0.882. 

Such  .a  break  in  the  curves  indicates  that  there  is  a discontinuity  in 
the  second  derivative  of  of  and  that  we  have  computed  a second  order 
transition  for  the  (100)  APB  with  classical  exponents. 

Another  way  of  describing  the  transition  is  in  terms  of  the  thickness 
of  the  (2,2)  layer  in  the  central  region  of  the  APB.  In  Fig.  12a  it  is 
tvn  atom  planes  thick  and  is  really  a disordered  (4)  layer.  It  becomes 


six  planes  thick  whan  KT/|w|  is  0,99  and  continues  to  increase  to  what 
should  be  infinite  thickness  at  the  ocnsolute  temperatures.  The  increases  in 
thickness  are  Monitored  in  Fig.  14  which  shows  two  curves  for  the  difference 
in  occupation  of  two  sublattices  on  a single  atom  plane  for  n=28  and  n=30 
(when  the  center  of  the  APB  occurs  between  n=30  and  n=31) . When  the 
difference  in  occupation  drops  to  zero  that  plane  has  joined  the  (2,2) 
layer.  The  n*30  curve  shows  the  disappearance  of  the  last  remnant  of  the 
(2,2)  layer  at  KT/|w|  = 0.882.  The  n=28  curve  shows  what  appears  to  be 
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an  abrupt  thicknening  near  KT/|w|  = 0.99.  There  is  a possibility  that  there 
are  a number  of  phase  transitions  associated  with  the  thickening  of  the 
(2,2)  layer.  It  was  not  possible  to  improve  the  precision  of  Fig.  13 
sufficiently  to  ascertain  whether  there  were  discontinuities  in  derivatives 
of  o'  corresponding  to  the  thickening  of  the  (2,2)  layer.  The  question  of 
these  additional  surface  phase  transition  at  the  higher  temperatures 
requires  further  study. 

Symmetry  requires  that  the  order  parameter  at  the  center  of  an  APB 
in  the  bcc  with  CsCl  order  go  to  zero.  No  such  symmetry  requirement 
exists  far  APB's  in  fee  with  Ll^  order,  and  APB's  below  the  0.882  transition 
were  ordered  at  thair  center.  Above  the  .882  transition  the  central 
two  planes  are  disordered  consistent  with  the  symmetry  of  an  APB  within 
a C2,2)  structure. 

The  structure  of  the  (11Q)  IPB  at  the  congruent  point  is  shown  in 
Fig.  15.  The  (110)  APB  is  not  shown  but  it  very  closely  resembles  two 
IPB's.  If  there  were  a C2,2)  layer  in  this  boundary  it  would  not  be 
apparent  in  Fig.  15  which  averages  compositions  on  (110)  planes  that  cut  at 
45°  to  the  layering  in  the  CuAu  structure.  We  did  not  investigate  whether 
or  not  a phase  transition  occurs  in  (110)  APB?s. 


Discussion 


Because  statistical  mechanical  calculations  of  foe  phase  diagrams 
have  only  recently  become  possible,  this  is  the  first  calculation  of 
equilibrium  properties  of  interfaces  within  and  between  phases  based 
on  fee.  The  first-order  character  of  the  phase  transitions  implies  that 
rwne  of  the  higher-order  critical  point  behavior  could  be  expected  for 
these  interfaces  at  the  disordering  transitions.  Indeed  the  surface  free 
energies  renain  finite  at  the  first-order  transition  point.  Nonetheless 
several  interesting  phenomena  were  uncovered.  The  (100)  APB's  undergo 
second  order  surface  transitions  in  which  the  internal  structure  of  the 
boundary  changes,  and  at  the  disordering  temperature  the  (h,k,0)  APB's 
contain  a thick  disordered  layer.  Indeed  the  APB  has  become  two  IPB's 
at  the  congruent  temperature  T,,  and  near  there  the  APB  is  perfectly  wet 
by  the  disordered  phase.  Modem  high  resolution  microscopy  should  be  able 
to  resolve  these  internal  structures. 

The  large  anisotropy  in  IPB  energy  expected  from  the  cubical  ly 
distorted  sphere  of  large  coherent  particles  was  not  calculated.  Indeed 
the  observation  that  snail  particles  are  more  spherical  is  consistent  with 
these  calculations  and  with  the  expectation  that  surface  energy  is  the 
dominant  determiner  of  stape  in  small  particle,  while  other  factors  such 
as  elastic  anisotropy  became  import  amt  for  large  particles.  APB's  were 
shown  in  Section  II  to  have  tetragonal  symmetry  and  to  be  surprisingly  close 
to  isotropic  in  the  zone  of  (001).  Because  of  our  neglect  of  second  neighbor 
interaction,  (001)  APB's  have  zero  energy.  This  leads  to  a very  high 
anisotropy  far  all  other  orientations  and  a prediction  that  all  (h,k,£)  APB's 
would  facet  into  (h,k,0)  and  (001).  Observations  on  APB's  indicate  both 
funded  and  facetted  portions,  and  such  domain  structures  should  be  reexamined 
to  establish  the  nature  of  the  difference. 


Me  next  turn  to  epplioability  of  our  results  to  reel  systems.  This 
hir^es  on  two  aspects  of  the  calculation,  the  model  vfdch  is  implicit  in 
tJse  expression  of  the  energy  in  terms  of  tetrahedral  clusters  in  equation (.5) 
and  the  tetrahedron  appraximationn  of  the  C.V.M.  It  is  important  to  distinguish 
%toat  ws  assume  in  the  model  for  the  energy*  which  is  couched  in  energies 
assigned  to  clusters*  with  the  statistical  problems  of  the  basic  clusters 
in  the  C.V.M.  In  this  calculation  we  used  energies  of  tetrahedra  for  the 
model  and  tetrahedra  for  the  C.V.M.  For  obvious  reasons  the  cluster  in  the 
C.V.M.  must  contain  (be  the  same  or  larger  than)  the  cluster  in  the  model, 
unless  one  adds  a Brogg-Williams  approximation  [20]. 

Mb  have  already  noted  that  for  fee  using  clusters  smaller  than  tetrahedra  in 
the  C.V.M.  gives  unrealistic  phase  diagrams,  while  using  tetrahedral 
clusters  in  a near-neighbor  pair  wise  energy  model  gives  a quite  realistic 
diagram  in  which  there  are  three  ordered  phases  which  however  have  three 
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congruent  points  / almost  the  same  temperature.  The  a and  g parameters 
represent  four-atom  near-neighbor  interaction  energies  which  create  differences 
in  congruent  point  temperatures.  The'  small  values  used  here  were  chosen  to 
fit  the  CuAu  diagram.  They  should  not  affect  the  qualitative  conclusion 
reached  here. 

Because  the  model  ignores  the  contribution  of  more  distant  neighbors  to 
the  interaction  energy  it  exaggerates  the  anisotropy  of  the  APB's . Oar 
finding  thatCf(OOl)  is  identically  zero  is  a direct  result  of  a near-neighbor 
model.  Introducing  higher  neighbor  interactions  would  require  an  increase 
in  the  size  of  the  cluster  in  the  C.V.M.  with  a large  increase  in  the 
computational  effort. 

The  accuracy  of  the  C.V.M.  is  usually  tested  by  having  seme  rigorous 
results  available.  For  bcc  APB's  such  results  existed  at  T=0  and  the  critical 
temperature,  and  it  was  possible  to  show  that  the  pair  approximation  gave 
quite  adequate  results.  Fbr  foe  the  results  merge  with  the  rigorous  predictions 
at  T=0. 


The  thermodynamic  self  consistency  checks  (Figs.  9 and  10)  test  neither 
the  validity  of  the  model  nor  the  approximation  in  the  C.V.M.  but  only  checks 
that  the  calculations  were  performed  consistently  with  both.  It  is  surprising 
how  many  calculations  on  surfaces  fail  this  check.  It  fails  e.g.  whenever 
there  is  an  artificial  constraint  on  the  thickness  of  the  boundary,  which 
prevents  it  fran  equilibrating  or  whenever  the  boundary  is  between  phases 
that  are  not  in  equilibrium. 

Corresponding  to  this  theoretical  check  is  the  experimented  use  of  the 
Gibbs  Adsorption  equation  as  a true  test  for  surface  equilibria.  An  extreme 
case  of  a nonequilibrium  surface  between  layers  of  water  and  sulfuric  acid 
gives  no  detectable  surface  tension.  Far  such  miscible  phases  one  would 
calculate  a negative  Cf  much  as  was  done  for  an  earlier  theory  of  IPB's  (23). 

Sane  results  of  our  calculation  are  completely  model  independent.  The 
large  difference  in  IPB  energy  on  either  side  of  the  congruent  point  and 
the  infinite  value  of  dO/dT  there  follows  directly  fran  the  Gibbs  adsorption 
equation.  It  implies  that  the  hysteresis  in  order-disorder  kinetics  due 
to  nucleation  will  be  quite  different  on  the  two  sides  of  any  congruent 
point.  For  the  same  reason,  coarsening  kinetics  of  ordered  precipitates 
driven  by  reduction  in  surface  area  would  be  quite  different  on  the  two 
sides  of  the  congruent  point. 
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TABLE  1 Definition  of  the  Pair  Variables 


Pair  of  Species  Probability  variables 


*pn  ” ^Qn 

Yn  (iJ) 

Vn  ’ kP(n^) 

Vppn  (l.k) 

*Pn  " kQ(n+l ) 

%.  <'•»> 

JQn  " fyn+l) 

V «-k> 

jQn  ' fyn+l ) 

VQQn 

TABLE  2 Definition  of 

the  Point  Variables 

Configuration  of  a Point  Probability  variables 


Figure  Captions 

An  antiphase  domain  boundary  separates  two  domains  of  the  same  ordered 
phase,  a)  (100)  projection  of  two  Li-  phases  showing  difference  between 
(010)  and  (001)  boundaries,  b)  (001)  projection  showing  that  (100)  and 
(010)  boundaries  are  quite  similar.  Excess  of  either  species  can  be 
accomodated  at  the  boundaries. 

The  convention  for  defining  the  four  sublattices  I,  II,  III  and  IV,  the 
planes  (...  n,  n+1,  ...),  the  layers  (P,Q) , and  the  tetrahedra. 

The  unit  cell  at  the  center  belong  to  one  domain  and  the  other  chree  belong 
to  another.  The  three- fold  axis  has  been  destroyed  and  two  of  the  three 
four-fold  axes  have  been  reduced  to  two-fold  axes.  Only  the  (001)  axis 
remins  a canton  four-fold  axis  for  both  domains. 

The  relationships  among  O',  the  ground  state  energy  of  (100)  APB's,  the 
composition  and  the  chemical  potential  is  represented  by  a single  curve. 

Its  projection  gives  the  relationships  among  any  two  of  them.  Thus,  while 
or  is  zero  except  at  stoichiometry  it  is  a continuous  function  of  u . 


The  Cu 
oonver 


Cu^Au 

ergerv 


Au  phase  diagram  showing  curves  of  constant  p 's. 
ence  at  n=Q.25. 


Note  the  non-uniform 


& of  QOO)  antiphase  boundaries  as  a function  of  temperature.  Curves 
are  either  for  constant  chemical  potential  u or  constant  composition  x. 

A curve  for  a (110)  APB  at  oonstant  composition  is  also  given. 

The  reduced  free  energy  O of  the  interphase  boundary  between  fee  and  the 
Ll-  structure  as  a function  of  reduced  temperature  for  (100)  and  (110).  The 
peak  in  O'  does  not  occur  at  the  peak  temperature . 

The  same  free  energy  as  Fig.  7 as  a function  of  reduced  chemical  potential 
instead  of  temperature. 

Thermodynamic  consistency  by  plotting  for  each  temperature  the  calculated 
value  of  O’  and  the  slope  calculated  from  equation  (18)  for  the  (100)  IPB. 
Canpare  with  Fig.  7. 

Consistency  with  the  Gibbs  adsorption  equation  for  APB's  is  demonstrated  by 
giving  for  each  temperature  the  calculated  values  of  Cf  and  its  temperature 
coefficient.  Compare  with  Fig.  6. 

The  calculated  profiles  of  the  (100)  APB  and  IPB  at  the  congruent  point  shows 
that  the  APB  has  become  two  IPB's.  The  composition  of  each  sublattice  is 
given  for  each  plane  parallel  to  the  interface.  Note  the  portion  where  the  two 
sublattices  in  each  plane  are  equally  occupied  (a  (2,2)  or  CuAu  layer). 

There  is  a structural  change  in  the  care  of  (100)  APB's  between  KT/|w|  = 0.90 
and  0.87. 


The  excess  entropy  and  adsorption  of  a (100)  APB  as  a function  of  temperature 
shows  that  the  transition  at  KT/|w|  = 0.882  is  second  order. 

The  ordering  of  atoms  on  the  two  sublattices  of  a given  (100)  atom  plane  as  a 
function  of  temperature  for  one  of  the  central  planes  (n=30)  and  the  third 
plane  from  the  center  (n=28). 
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Figure  la. 
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Figure  lb. 
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Recent  calculations  of  a two-dimensional  lattice  gas  model  has 

i o 

led  to  vapor,  liquid  and  crystalline  phases.  * In  addition  to 
permitting  the  phase  transitions  to  be  studied,  such  a model  can  be 
used  in  studying  linear  interfaces  between  all  pairs  of  phases  as 
well  as  the  grain  boundary  between  two  crystals  of  the  same  solid 
phase  differing  in  orientation  (Figure  1).  We  report  finding  a gradual 
but  well  defined  transition  in  the  grain  boundary  structure  well 
below  the  melting  point  of  the  crystal;  as  the  transition  region  is 
passed,  the  grain  boundary  structure  changes  from  its  low  temperature 
configuration  into  such  a structure  that  there  is  a layer  along  the 
grain  boundary  whose  thickness  increases  to  infinity  as  the  logarithm 
of  the  undercooling  and  that  the  properties  of  the  grain  boundary  - 
approach  that  of  two  solid-liquid  interfaces  separated  by  a liquid. 

As  the  melting  point  is  approached,  the  grain  boundary  shows  a 
singular  behavior  in  that  the  excess  entropy  and  the  excess  specific 
heat  due  to  the  grain  boundary  (of  a unit  length)  becomes  infinite. 

Our  model  is  a two-dimensional  square  lattice  gas  with  inter- 
action potentials  chosen  to  match  a prior  computer  simulation1.  We 
used  the  cluster  variation  method  (CVM)  and  natural  iteration  (NI) 
to  calculate  the  equation  of  state  and  the  phase  diagram.  Figure  2 
indicates  that  the  phase  diagram  so  calculated  matches  that  obtained 
by  the  computer  study.  The  properties  of  interphase  interfaces  and 
the  grain  boundary  shown  In  Fig.  1 were  calculated  along  the  appro- 
priate two-phase  coexistence  of  Fig.  2 and  for  the  grain  boundary 


between  two  solid  state  domains  of  different  orientations  as  In 

Fig.  1;  the  CVM  and  NI  were  used  as  In  our  previous  studies  of  antl- 

5-  7 

phase  domains  boundaries  In  ordered  structures.  If  we  consider 

crystallization  model  to  be  an  ordering  of  holes  (V)  and  atoms  (A) 

into  a V^A  structure,  the  grain  boundary  Is  Indeed  a boundary  between 

two  ordered  domains.  It  is  also  a model  of  domains  In  an  adsorbate 

layer  which  can  crystallize  with  a rotation  with  respect  to  the 

8 9 

underlying  substrate.  * 

The  particular  grain  boundary  In  Fig.  1 Is  a symnetric  tilt 
boundary.  With  this  much  tilt  a I B 5 coincidence  lattice  exists 
in  which  one  In  five  atoms  (e.g.,  P and  P')  occupy  a site  that  could 
belong  to  either  crystal  structure.  The  underlying  lattice  Is  known 
as  the  DSC  lattice.10  If  we  let  a be  the  lattice  constant  of  the 
DSC  lattice,  the  lattice  constant  of  the  crystal  (as  marked  by  white 
or  black  circles  In  Fig.  1)  Is  Js  a. 

For  the  basic  cluster  In  CVM,  we  used  the  nine-site  cluster 
shown  at  the  bottom  of  Fig.  1.  Pair  of  atoms  were  not  allowed  closer 
than  AB,  the  energy  of  pairs  at  distances  AB  and  AC  was  assumed 
attractive,  being  respectively  -1.2e  and  There  is  no  Inter- 

action between  atoms  further  apart  than  the  distance  AC.  We  calculate 
the  equilibrium  state  of  the  entire  system  including  the  boundary,  and 
obtain  the  probabilities  of  encountering  each  of  twenty  possible 
arrangements  of  atoms  on  the  nine-site  clusters  centered  on  each  site 
in  a strip  of  width  PP*  along  the  grain  boundary  and  extending  from 


j - 1 to  j - 81  (sometimes  61  or  41)  Into  each  grain.  From  these 
cluster  probabilities  we  compute  the  excess  (compared  to  an  equili- 
brated single  crystal)  In  the  grand  potential  for  this  strip 
o = E - TS  - yNa  (1) 

- where  E is  the  excess  potential  energy,  T Is  the  absolute  temperature, 

S Is  the  excess  entropy,  u is  the  chemical  potential  and  N is  the 

d 

excess  number  of  atoms  in  the  strip  of  width  PP'.  We  will  define  the 
unit  length  along  the  boundary  to  be  the  distance  PP'  = S\/2  a where 
a is  the  lattice  constant  of  the  underlying  DSC  lattice.  Thermodynamic 
self-consistency  requires  that11 

do  = -SdT  - N dp  (2) 

which  is  a useful  test  for  equilibrium  and  forms  the  basis  for  our 
extrapolation  (v.i.). 

Fig.  3 shows  o calculated  as  a function  of  temperature  for 
p/r  = -1.5.  The  low  temperature  behavior  is  readily  understood  by 
examining  Fig.  1.  The  pair  QQ'  are  forbidden  and  either  site  must 
be  empty.  For  the  unit  length  of  the  boundary  this  leads  to 
S = k ln2,  Ng  = -1  and  E = 5^  - at  T ■ 0.  Hence  for  low  o 
o =5e  + p-  kT  ln2  (3) 

which  agrees  with  the  curve  for  kT/e<0.3. 

The  density  profiles  perpendicular  to  the  boundary  are  shown  in 
Fig.  4 for  various  temperatures.  The  low  temperature  W shape  near 
the  center  of  the  boundary  is  consistent  with  the  expectation  that 
the  layer  on  either  side  of  center  is  half  occupied.  As  the  melting 


* 


J 


L- 


point  Is  approached  a low  density  layer  Is  formed  near  the  center  of 
the  boundary  and  the  thickness  of  the  layer  approaches  our  computer 


capacity.  The  excess  quantities  S and  -N  tended  to  Increase  without 

d 

limit  as  did  -3c/ 3T  as  the  melting  point  Is  approached. 

Figure  5 plots  S vs.  -log(T  -T).  For  y/e  - -1.5,  T Is  0.71629  e/k. 

m m 

It  Is  observed  that  the  curve  In  Figure  5 Is  made  of  two  distinctly 
different  portions.  For  low  temperatures  (kT/e<0.20),  we  see 
S - k ln2  (4) 

for  high  temperatures  (kT/e>0.45),  S Is  linear  In  log(Tm-T)  as 
S - -3.681  - 4.1520  ln(Tm-T)  (5) 

Making  use  of  these  relations  for  S In  Figure  5,  we  can  Integrate 

/sdT  to  obtain  the  estimate  of  o(T)  for  the  high  temperature 

region: 

oHJ  - 1.523  - (Tm-T)  -0.471  + 4.152  ln(Tm-T)  (6) 

This  Is  shown  as  a solid  curve  In  Fig.  3. 

The  good  agreement  between  the  solid  curve  crHT  and  Individually 
computed  o values  Indicated  by  dots  shown  In  Fig.  3 serves  as  a test 
of  Eg.  (2),  the  self-consistency  of  the  formulation,  but  It  also 
permits  extrapolation  of  o to  the  melting  temperature  where  its  value 
Is  1.523. 

We  separately  calculated  the  solid-melt  Interfacial  properties 
for  the  orientation  corresponding  to  that  of  our  grain  boundary.  The 
value  of  Its  reduced  o was  found  to  be  0.763,  while  Its  density 
profile  closely  matched  that  of  either  half  of  the  grain  boundary. 


We  concluded  not  only  that  the  grain  boundary  has  become  a melted 
layer  at  the  melting  point,  but  that  it  behaves  as  If  it  is  coated 
with  a melted  layer  for  a considerable  temperature  interval  below 
the  melting  point.  At  some  temperature  near  0.35c/k  there  is  a gradual 
transition  in  structure  from  the  low  temperature  profile  to  a struc- 
ture which  with  Increasing  temperature  increasingly  tends  to  resemble 
a melted  layer  and  which  exhibits  a singularity  as  Tm  is  approached. 
Since  the  grain  boundary  is  completely  wet  at  T , it  is  legitimate  to 
call  the  gradual  transition  near  T = 0.35e/k  the  wetting  transition. 

The  possibility  that  a grain  boundary  would  have  a liquid  layer 

12 

at  the  melting  point  was  first  discussed  by  Gibbs.  Whenever  a 

for  the  grain  boundary  has  a tendency  to  exceed  twice  the  o for  the 

solid-liquid  interface  the  transition  we  found  is  expected  to  occur. 

13 

Smith  attempted  to  formulate  the  temperature  behavior  of  this 

transition  by  comparing  two  models  of  the  grain  boundary.  The  solid 

model  was  assumed  to  possess  a known  value  of  a.  The  value  of  a 

for  the  melted  layer  model  was  assumed  to  consist  if  three  terms 

c^L  * 2c^l  + XAS(Tm-T)  + E(A)  (7) 

where  2a^  is  the  contribution  of  the  two  solid-melt  interfaces, 

AAS(T  -T)  is  the  contribution  of  a melted  layer  of  thickness  A,  AS 
m 

is  the  entropy  of  fusion,  and  E( A)  is  an  unknown  repulsive  energy. 

If  2o$l  was  less  than  o,  the  melted  layer  took  over.  Minimizing 

with  respect  to  A at  a fixed  temperature  we  obtain  a relation  for  A 
DE/dA  + AS(Tm-T)  = 0 


(8) 


When  this  relation  holds,  and  when  °Sl  and  AS  are  assumed  Independent 
of  temperature,  we  can  further  obtain 

doML/dT  * 'XAS 

Since  this  derivative  Is  -S,  our  finding  In  Fig.  5 and  Eq.  (5)  Implies 
that  X linearly  depends  on  -ln(Tm-T),  and  integration  of  Eq.  (8)  leads 
to  an  expression: 


E(X)  * E(0)  + C^xpt-CgX)  (10) 

This  Indicates  the  reasonable  nature  of  the  repulsive  term  in  Smith's 
formulation  (7),  When  we  use  these  relations,  we  can  write  the 
tenperature  dependence  of  In  (7)  in  the  form 


°ML=  2°Sl  ' <VT) 
which  has  the  name  (T  - 

W 


c3  * C4  1»< T„-T) 


(11) 


) -dependence  as  (6). 

In  the  present  paper,  we  used  a two-dimensional  lattice  gas- 
liquid-solid  model,  and  thus  the  grain  boundary  is  essentially 
one-dimensional.  This  one-dimensional  behavior  is  consistent  with  the 
nature  of  the  gradual  transition  from  the  low- temperature  behavior 
to  the  high-temperature  behavior  demonstrated  clearly  In  Fig.  5.  If 
we  work  on  a three-dimensional  system  with  a two-dimensional  grain 
boundary.  It  is  expected  that  the  nature  of  the  low-  to  high- temperature 
behavior  may  be  a sharper  one;  for  example,  a second-order  phase  transition. 

A melting  transition  in  real  metals  has  been  observed14  as  a function 
of  orientation  difference  and  explained  in  terms  of  a dislocation  model 
with  liquid  cores. ^ The  temperature  dependence  of  such  a model  might 


show  several  transitions. 
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Figure  1.  A two-dimensional  model  of  a grain  boundary. 

Position  of  atoms  in  the  crystalline  state 
are  shown  by  circles.  This  illustration 
is  the  unrelaxed  state,  and  thus  is 
schematic . 
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Figure  4(a).  Density  profile  across  the  grain  boundary.  The  center 
of  the  boundary  is  chosen  at  the  41st  lattice  plane. 
The  profile  is  repeated  on  the  right  of  the  41st  plane 
as  a mirror  image,  p/e  ■ -1.5  and  starts  from  a very 
low  temperature. 
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Figure  4(b).  Density  profile  across  the  grain  boundary.  The  center 
of  the  boundary  is  chosen  at  the  41st  lattice  plane. 
The  profile  is  repeated  on  the  right  of  the  41st  plane 
as  a mirror  Image,  y/e  ■ -2.2  and  shows  how  the 
profile  changes  near  the  melting  temperature  (Tm  - 
0.5031638).  The  curve  for  Tm  is  for  an  intermediate 
stage  of  convergence. 


